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Abstract 

The  objective  of  the  present  research  is  to  compute  the  variation  of  tur¬ 
bulent  Prandtl  number  across  the  viscous  sublayer  for  fluids  of  different  molecular 
Prandtl  number  and  for  flows  with  adverse  and  favorable  pressure  gradients.  This 
is  accomplished  using  a  time-dependent  compressible  Navier-Stokes  code  together 
with  a  coherent  structure  model  of  the  boundary  conditions.  Existing  theories  for 
turbulent  Prandtl  number  Prr  near  a  wall  are  shown  to  differ  greatly,  while  exper¬ 
imental  measurements  are  scarce  and  exhibit  a  large  uncertainty  band. 

A  variety  of  boundary  conditions  for  the  coherent  structure  model  were 
explored  and  are  briefly  described.  Numerical  results  are  presented  for  the  set  of 
boundary  conditions  that  yielded  the  most  realistic  results  for  those  turbulence 
quantities  which  could  be  compared  with  experiment. 

Computational  results  for  Prx  have  been  obtained  for  the  cases  of  varying 
molecular  Prandtl  number  from  0.72  to  6.0;  for  pressure  gradients  in  wall  variables 
of  ^0.01  (mild  adverse), zero,  and  -0.01  (mild  favorable);  and  for  two  different  wall 
temperature  boundary  conditions.  These  results  are  compared  to  various  theories 
and  experiments.  It  is  found  that  the  dependence  of  PrT  on  pressure  gradient  is 
negligible;  but  that  the  dependence  on  molecular  Prandtl  number  is  sizable  near 
the  wall  [y*  <  3).  Turbulent  Prandtl  number  is  found  to  be  essentially  constant 
throughout  most  of  the  viscous  sublayer  for  the  Prandtl  number  and  pressure  gra¬ 
dient  variations  investigated.  It  is  shown  that  the  assumption  of  constant  Ptj  (0.9) 
across  the  viscous  sublayer  is  adequate  for  practical  heat  transfer  computations. 
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1.  Introduction  and  Objectives 

Many  engineering  problems  today  are  concerned  with  some  aspect  of  tur¬ 
bulent  flow,  and  are  complicated  by  the  complex  nature  of  turbulence.  For  most 
calculations,  some  type  of  modeling  is  used  since  it  would  take  huge  amounts  of 
computer  time  to  capture  all  of  the  turbulent  dynamics  numerically.  Many  calcu¬ 
lations  use  Reynolds'  averaging.  Once  this  is  done  however  the  set  of  conservation 
equations  for  mass. momentum,  and  energy  is  not  closed.  This  closure  problem  is 
solved  by  expressing  the  added  unknowns  in  terms  of  the  averaged  flow  variables. 

Many  of  the  terms  which  need  to  be  modeled  for  turbulent  flow  are  dominant 
in  the  viscous  sublayer  region  very  close  to  the  wall.  Turbulence  kinetic  energy 
production  and  dissipation  have  their  peak  values  in  the  viscous  sublayer  at  about 
a  y~  10  (Hinze  1975).  Also  70T  of  the  production  of  turbulence  occurs  in  the 
viscous  sublayer  (Kim  et.  al.  1971).  In  the  case  of  turbulent  heat  transfer,  the  value 
of  turbulent  Prandtl  number  in  the  viscous  sublayer  has  been  found  to  be  a  very 
important  parameter  in  determining  the  rate  of  heat  transfer  and  the  temperature 
profile.  The  turbulent  Prandtl  number  appears  when  the  energy  equation  is  time- 
averaged.  and  is  defined  as 


Prj  = 


Difficulties  are  encountered  in  measuring  the  turbulent  Prandtl  number  near 
wall  in  the  viscous  sublayer  (hereafter  called  YSL).  This  is  because  experimental 
measurements  are  technically  difficult  to  obtain  due  to  the  extreme  thinness  of  the 
viscous  sublayer.  For  a  high  Reynolds  number.  Re i  =  2.0x10*'.  wing  cord  length 
L  -  10m  and  free  stream  velocity  Cx  ~  300m  s  in  air.  the  viscous  sublayer  thick¬ 
ness  dui  -  .03mm.  This  is  approximately  .05c_c  of  the  boundary  layer  thickness  t. 
and  is  typical  for  a  boundary  layer  on  a  jet  transport  wing.  For  a  low  Reynolds 
number.  Rei  ~  1  0x10'  .  a  fiat  plate  length  L  -  1.0m.  and  free  stream  velocity 
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l'x  —  20.0m  s  in  air.  the  viscous  sublayer  thickness  blf{  -  ,7mm.  This  is  approx¬ 
imately  5^  of  the  boundary  layer  thickness  6.  This  is  typical  for  an  experimental 
channel  flow  with  air.  For  the  higher  Re.  it  is  clear  that  the  YSL  is  very  thin.  The 
thickness  of  the  viscous  sublayer  can  be  enhanced  with  low  Re  and  fluids  of  high 
viscosity.  In  wall  units,  the  thickness  of  the  YSL  is  about  y~  -  40.  Wall  units 
are  the  conventionally  used  nondimensionalization  in  turbulence  analysis.  The  YSL 
designates  the  region  from  the  wall  to  the  beginning  of  the  log  region  which  is  at 
about  y~  ^  30  -  50.  For  this  study,  the  edge  of  the  YSL  will  be  at  40.  Because  the 
YSL  is  so  thin,  near  wall  measurements  of  tZr  and  vO  correlations  are  unreadable 
due  to  probe  interference  effects.  This  is  because  the  long-dimension  of  the  probe 
has  to  be  normal  to  the  wall  in  order  to  measure  the  v-velocity  component.  The 
measurements  of  the  normal  gradients  of  u  and  9  are  difficult  because  within  the 
thin  YSL.  these  gradients  are  very  steep.  Hence,  measurements  of  Prr  near  a  wall 
are  highly  uncertain. 

Some  experimentalists  use  alternate  formulations  for  the  turbulent  Prandtl 
number  which  involve  additional  assumptions.  These  still  rely  on  having  accurate 
values  of  the  velocity  and  temperature  gradients  near  the  wall  so  their  reliabilty 
can  still  be  questioned.  An  example  of  this  is  the  expression  used  by  L.C.  Thomas 
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36  "  Pr 
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(1.2o) 


which  assumes 


The  extreme  experimental  difficulties  in  measuring  Pr j  make  all  measure¬ 
ments  of  PrT  in  the  YSL  unreliable.  The  experimental  uncertainty  envelope  for  air 
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with  Fr  =  .72  is  shown  in  figure  1  (Simpson  et.  al  1973).  This  shows  that  even 
for  air  as  the  wall  is  approached  the  scatter  of  the  measured  value  of  the  turbulent 
Prandtl  number  is  large.  Launder  (1978)  states  that  for  Fr  •  1.  the  scatter  of 
experimental  data  is  so  large  even  awa\  from  the  wall  that  the  existing  data  is  hard 
to  evaluate  except  to  sa\  there  is  much  scatter.  Experimental  data  for  Fr  1  and 
of  non-zero  pressure  gradients  for  even  casual  perusal  is  scarce.  Figure  2  shows  two 
data  sets  for  the  Fr with  adverse  pressure  gradient.  Thomas’  data  (19821  i*-  from 
his  computational  model.  Blackwell’s  data  (1972)  is  experimental.  These  two  data 
'et^  show  contrasting  behavior  of  the  Fr  j  in  an  adverse  pressure  gradient 

Many  analytic  expressions  have  been  developed  for  the  turbulent  Prandtl 
number,  as  illustrated  in  the  survey  of  Reynolds  (1975).  The  initial  assumptions  of 
these  various  models  of  the  turbulent  Prandtl  number  differ  greatly.  A  basic  list  of 
the  models  would  be  the  following: 

(1)  Various  mixing  length  models  like  the  van  Driest  temperature  wall 
damping  function  from  Stokes  flow  analogy 

(2)  Turbulent  heat  transfer  models  like  Jenkin’s  or  Diessler’s  analysis 
with  heat  transfer  to  and  from  moving  spherical  edd\ 

(31  Formal  analysis  based  on  exact  equations  like  Geshev’s  Green's 
function  analysis  on  the  linearized  equations 

The  validity  of  these  models  is  hard  to  estimate  since  the  experimental  data 
varies  tremendously.  Various  existing  theories  are  shown  in  figures  3-6  for  different 
molecular  Prandtl  numbers  between  0.72  and  1000.  As  with  the  experimental  data, 
the  turbulent  Prandtl  number  exhibits  wide  scatter  as  the  wall  is  approached.  The 
scatter  of  values  of  Ptj  is  not  unexpected  since  the  experimental  measurements  of 
PrT  can  give  no  trends  of  its  behavior  in  the  V5L  upon  which  a  reliable  theory  can 
be  based.  For  Pr  -  .72.  the  scatter  is  not  as  wide  since  most  of  the  expressions 
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derived  for  the  turbulent  Prandtl  number  include  at  least  one  floating  coefficient 
which  was  adjusted  to  give  a  Pr j  of  order  1  in  the  outer  part  of  the  VSL. 

The  turbulent  Prandtl  number  is  very  important  in  the  determination  of 
heat  transfer:  the  existing  state  of  knowledge,  however,  is  incomplete  and  confusing. 
The  present  theories  and  experimental  data  sets  show  both  strong  qualitative  and 
quantitative  differences  which  are  impossible  to  collapse  onto  any  unified  curve. 
Several  general  conclusions  have  been  advanced  from  research  on  the  turbulent 
Prandtl  number. (Reynolds  1975.  Jischa  1978) 

(1 )  Ft  :  1  then  Prj  <  1 

(2)  Pr  <r  1  then  Pr j  >  1 

(3)  <3-^—  >  0  then  Prj  decreases  from  its  =  0  value 

ai~  ai~ 

These  conclusions  only  apply  to  the  region  away  from  the  wall.  There  are  no  general 
conclusions  for  the  behavior  near  a  wall. 

This  present  state  of  great  uncertainty  about  the  near-wall  behavior  of  the 
turbulent  Prandtl  number  gives  an  obvious  motivation  to  u  idertake  a  computational 
investigation.  The  objectives  of  the  present  research  are  to  study  the  behavior  of 
the  turbulent  Prandtl  number  as  a  function  of  various  parameters.  Specifics  of 
the  study  are  to  investigate:  (1)  the  dependence  of  the  PrT  on  molecular  Prandtl 
number  between  0.72  ■  f‘r  •  6.0:  (2)  the  effect  of  a  mild  mean  streamwise  pressure 
gradient  both  adverse  and  favorable  on  Prj\  and  (3)  the  effect  on  Prj  of  using 
different  temperature  wall  boundary  conditions.  Also  investigated  is  the  limiting 
behavior  of  Prj  as  the  wall  is  approached. 

The  method  chosen  for  this  study  is  that  of  a  time  dependent  three  dimen¬ 
sional  Navier  Stokes  computation  with  coherent  structure  modeling.  This  approach 
was  developed  by  Chapman  and  Kuhn  (1981)  and  Hatziavramidis  and  Hanratty 
/ 1979)  for  modeling  viscous  sublayer  turbulence  in  incompressible  flow  without  heat 
transfer  .  During  the  past  two  decades  a  great  deal  of  new  experimental  information 
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has  been  assembled  on  the  physics  of  organized  eddy  structures  in  turbulent  flow 
especially  within  the  sublayer  (Cantwell  1981).  Yet  it  has  not  been  possible  thus  far 
to  incorporate  this  body  of  physical  information  within  the  framework  of  Reynolds- 
averaged  turbulence  modeling.  A  fundamental  limitation  of  the  Reynolds-averaged 
approach  is  that  time  averaging  is  done  to  the  dynamic  equations  of  motion.  In  this 
initial  mathematical  step  important  physical  aspects  of  organized  eddy  structures, 
such  as  phase  relationships  and  coherent  structure  dynamics  are  obliterated  irre¬ 
versibly.  Consequently  some  totally  different  approach  is  required  if  the  observed 
physics  of  coherent  eddy  motions  are  to  be  incorporated  within  the  framework  of  a 
t urbulence  model. 

Coherent  structure  modeling  fits  this  requirement.  This  method  directly 
models  the  essential  organized  eddy  structures  observed  in  experiments.  The  prin¬ 
ciple  steps  in  coherent  structure  modeling  are  the  following: 

(1)  model  velocity  and  temperature  boundary  conditions  at  the  edge 
of  the  viscous  sublayer 
(2 1  compute  time-dependent  dynamics 

(3)  time  average  computed  results 

Thus  time  averaging  is  the  last  operation  performed  on  the  computed  dynamics, 
rather  than  the  first  operation  performed  on  the  dynamic  equations.  The  study  by 
Chapman  and  Kuhn  (1981)  led  to  surprising  success  in  modeling  some  important 
characteristics  of  viscous  sublayer  turbulence. 

Both  the  motivations  and  the  payoff  for  the  present  research  are  quite  dif¬ 
ferent  from  those  in  conventional  turbulence  modeling.  In  addition  to  the  specific 
motivations  stated  for  the  present  research,  there  are  other  important  motivations 
that  could  affect  different  areas  of  turbulence  research.  It  is  appropriate,  therefore, 
to  outline  some  of  these. 


Three  other  motivations  for  developing  a  realistic  coherent  structure  model 
of  viscous  sublayer  turbulence  are  the  following: 

(1)  To  provide  a  basis  for  strengthening  a  major  ueak  link  m  present 
Reynolds  average  closure  schemes  by  providing  a  method  of  com¬ 
puting  the  production. destruction,  and  redistribution  terms  in  the 
exact  transport  equation  for  dissipation  or  turbulent  heat  flux  c\  hich 
then  can  be  compared  against  the  corresponding  modeled  terms. 

(2)  To  provide  a  simple  test  flow  against  which  subgrid  models  of  tur¬ 
bulence  in  large  eddy  simulations  can  be  tested.  A  realistic  coherent 
structure  model  of  viscous  sublayer  turbulence  could  be  used  to  test 
subgrid-scale  models  for  large  eddy  simulations  of  turbulent  flow  in 
the  region  adjacent  to  a  wall. 

(3)  To  provide  a  guide  for  modeling  the  lower  boundary  conditions  for 
the  outer  turbulent  region  in  large  eddy  simulations  of  boundary 
layer  flow.  If  a  realistic  boundary  condition  of  this  type  could  be 
developed,  the  entire  viscous  sublayer  could  be  modeled  rather  than 
directly  computed  which  would  greatly  help  large  eddy  simulations. 

In  summary,  with  these  motivations  for  developing  coherent-structure  YSL  turbu¬ 
lence  models,  it  is  clear  that  the  potential  payoffs  may  be  of  much  broader  scope 
than  the  specific  objectives  of  the  present  research. 


2.  Effect  of  PrT  on  Heat  Transfer 
The  turbulent  Prandtl  number  is  an  important  parameter  which  appears 

in  the  turbulent  energy  equation  when  the  time-average  is  taken.  Thus.  Prj  is  a 
parameter  which  must  be  modeled  in  computations  of  turbulent  heat  transfer  from 
the  turbulent  boundary  layer  equations. 

In  the  case  of  turbulent  heat  tranfer,  the  turblulent  heat  flux  is  rewritten 
using  a  Boussinesq  approximation.  This  still  involves  an  unknown  parameter,  the 
thermal  eddy  diffusivity,  tq. 

qT  -  -pcpv$  (2.1) 

—  , 

In  order  to  simplify  the  solution  a  little  more,  Reynolds  assumed  that  the  mecha¬ 
nisms  of  energy  and  momemtum  transfer  were  similar,  and  then  applied  dimensional 
analysis.  W  hen  this  is  done,  the  turbulent  Prandtl  number  is  related  to  the  thermal 
eddy  diffusivitv  in  a  manner  analogous  to  the  molecular  Prandtl  number.  Both  Pr 
and  PrT  are  ratios  of  measures  of  the  diffusivity  of  momentum  to  temperature. 


PrT  = 


<?r  =  -pc? 


tm  d6~ 
PrT  dy~ 


Conventionally,  the  eddy  viscosity  em  is  obtained  from  an  empirical  equation  or  a 
hydrodynamic  balance  equation.  The  turbulent  Prandtl  number  is  also  obtained 
from  an  empirical  equation.  Finally,  both  these  values  are  then  used  in  the  turbulent 
energy  equation  which  can  now  be  solved.  The  turbulent  Prandtl  number  is  used 
extensively  in  many  computations;  however,  its  value  for  fluids  of  varying  Pr  and 
its  behavior  under  different  flow  conditions  is  not  well  understood. 


Many  investigators  have  conducted  studies  on  the  turbulent  Prandtl  number 
with  the  main  objective  of  estimating  its  value  for  many  different  molecular  Prandtl 


numbers  and  types  of  flows.  From  Reynolds’  original  analysis,  the  turbulent  Prandtl 
number  is  a  constant  equal  to  1.  This  is  called  Reynolds  analogy.  The  effect  of  the 
presence  of  the  wall  was  known  to  be  an  added  difficulty.  Near  the  wall  viscous 
and  molecular  processes  were  thought  to  be  dominant  and  the  turbulent  Prandtl 
possibly  to  be  an  unimportant  parameter.  Some  experiments  (Johnson  1959)  have 
suggested  that  the  turbulent  Prandtl  number  may  not  be  constant  across  the  viscous 
sublayer.  Once  this  wras  observed,  several  researchers  investigated  the  effect  of  a 
variable  turbulent  Prandtl  number  on  the  flow  solution.  Rotta,  von  Driest,  and 
Reichardt  (Schlichting  1979)  all  concluded  that  the  value  of  the  turbulent  Prandtl 
number  near  wall  in  the  viscous  sublayer  is  an  important  parameter  in  determining 
the  rate  of  heat  transfer  and  the  temperature  profile. 

A  simple  computation  serves  to  show  this.  This  computation  is  similar 
to  that  in  Rotta’s  paper  (1964)  reformulated  to  emphasize  the  importance  of  the 
turbulent  Prandtl  number.  An  equation  of  temperature  as  a  function  of  y  +  .  distance 
from  the  wall,  Pr,  molecular  Prandtl  number,  Pr^,  turbulent  Prandtl  number,  and 
cm.  eddy  viscosity  can  be  written  as  follows  (Kays  and  Crawford  1980). 


-f -r 


dy' 


(2.5) 


Pr  v 

The  van  Driest  formulas  for  eddy  viscosity  and  the  mean  streamwise  velocity  gra¬ 
dient  are  (Hinze  1975) 


r, 

1  -  exp 


-y~  -  2  du  + 
A  J  dy+ 


(2.6) 


du~ 

dy'r 


V 


'  1  -4 «2y 1  -  exp 


-',2 


(2.7) 


where  A  is  a  fitted  parameter  which  has  a  value  of  26.  Then  for  a  given  Pr.  a 
two  layer  step  model  for  the  Pr-j-  is  used  to  assess  the  effects  on  heat  transfer.  The 
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interface  position  y,  between  the  two  steps  is  also  an  additional  variable.  The  results 
for  Pr  -  .72  (air)  and  Pr  -  6.0  (water)  are  shown  graphically  in  figures  7  and  8.  The 
calculation  was  first  made  for  a  base  value  of  the  Pr?  =  1.0.  Then  this  was  increased 
and  decreased  by  a  factor  of  two  in  the  inner  step  region  (y  <  y,).  The  interface 
position  y,  was  moved  out  until  the  approximate  heat  transfer  coefficients  differed 
by  approximately  20%  at  the  end  of  the  integration.  An  approximate  heat  transfer 
coefficient  was  computed  using  the  temperature  at  the  end  of  the  integration.  These 
values  were  used  to  get  an  approximate  change  in  the  heat  transfer  coefficient  A h 
in  per  cent.  Therefore  within  the  value  of  y,  shown,  the  value  of  the  Prj  does  not 
effect  the  heat  transfer  solution  significantly  (<  20%).  For  Pr  =  .72,  this  value  is 
around  y~  =  13.7,  and  for  Pr  =  6.0,  the  value  is  about  y+  =  7.5.  It  can  be  seen 
that  the  value  of  the  Ptj  is  important  even  in  regions  where  molecular  diffusion 
should  be  dominant.  This  region  moves  closer  to  the  wall  as  the  Pr  increases  which 
is  expected  since  the  thermal  layer  is  getting  thinner. 

Therefore,  if  quantitatively  accurate  heat  tranfer  computations  are  desired 
then  accurate  values  of  the  turbulent  Prandtl  number  in  the  viscous  sublayer  are 
needed  as  an  input.  Since  the  value  of  turbulent  Prandtl  number  in  the  VSL  in 
turbulent  heat  transfer  computations  is  important,  and  since  experiment  and  theory- 
have  not  been  able  to  define  the  variation  of  the  Prj  w-ith  molecular  Pr  and  mean 
streamwise  pressure  gradient  (as  previously  shown), there  is  a  need  for  a  Navier- 
Stokes  computation  of  the  PrT  distribution  in  the  viscous  sublayer. 


3.  Formulation  of  Boundary  Conditions 
for  the  Viscous  Sublayer 

The  study  by  Chapman  and  Kuhn  (1981)  showed  promising  viscous  sublayer 
turbulence  results  for  flow  with  zero  heat  transfer.  In  the  present  study,  a  similar 
relatively  simple  coherent-structure  model  is  developed  for  the  case  of  flow  with 
heat  transfer.  Thus  the  present  research  provides  an  extension  of  their  work  to  flow 
with  heat  transfer. 

The  primary  element  in  formulating  the  computational  model  involves  the 
construction  of  time  and  space  dependent  boundary  conditions  that  are  as  realistic 
as  possible  for  the  fluctuating  temperature  6  and  the  fluctuating  velocity  compo¬ 
nents  u.x-.u-  at  the  outer  edge  of  the  viscous  sublayer.  These  boundary  conditions 
are  constructed  using  various  experimental  observations  as  a  guide  to  model  the 
coherent  structures  present. 

The  boundary  conditions  attempt  to  model  the  bursting  process  revealed 
by  experimental  observations  of  flow  in  the  viscous  sublayer.  The  bursting  process 
comprises  two  events,  ejection  sweep  events  and  sweep  ejection  events  (Torino  and 
Brodke\  1969.  Johansson  and  Alfredsson  1982).  These  events  are  associated  with 
the  highly  elongate  streamwise  streaks  present  in  the  viscous  sublayer.  The  ejec¬ 
tion  event  is  associated  with  fluid  being  lifted  away  from  the  wall  from  a  low  speed 
streak  which  is  a  region  of  slower  moving  streamwise  fluid.  A  characteristic  of  the 
ejection  process  is  fluid  being  lifted  away  from  the  wall  and  creating  an  inflectional 
profile  of  the  streamwise  velocity.  When  fluid  moves  towards  the  wall,  high  speed 
streaks  which  are  regions  of  faster  moving  streamwise  fluid  occur.  This  is  called 
a  sweep  event.  Ensemble  conditional  averaging  shows  evidence  of  contra-rotating 
vortical  motion  associated  with  these  high  and  low  speed  streaks  (Blackwelder  and 
Eckelmann  1979).  Iritani  (  1981  )  has  shown,  using  liquid  crystal  flow  visualization 
and  the  comparison  of  temperature  correlations  to  velocity  correlations,  that  the 
temperature  field  under  slight  heating  conditions  also  has  evidence  of  a  vortical 


structure.  This  is  logical  since  temperature  should  be  acting  as  a  passive  con¬ 
taminant.  A  computational  study  by  Kim  (1984)  has  shown  that  the  the  vortical 
structure  does  not  have  any  “long  legs"  which  lie  along  the  wall,  but  instead  these 
legs  have  characteristic  lengths  much  shorter  than  the  typical  length  scale  of  the 
streaks.  These  vortex  legs  lift  off  immediately  forming  the  horeshoe  vortex  structure 
at  an  angle  of  45  degrees.  This  is  in  constrast  to  the  model  suggested  bv  Blackwelder 
and  Eckelmann  of  highly  elongate  contra-rotating  vortical  pairs  associated  with  the 
streaks  with  characteristic  lengths  of  the  order  of  1000  wall  units.  (  see  figure  9 
).  Moser  (19841  has  indicated,  also  from  a  computational  study.  that  the  vortical 
structure  does  not  occur  as  a  contra-rotating  pair  in  a  time-instantaneous  frame  but 
only  in  ensemble-averaging.  In  developing  the  present  YSL  model,  a  simple  ideal¬ 
ized  picture  of  the  time-instantaneous  contra-rotating  pairs  in  the  YSL  is  assumed. 
Regions  where  the  vortices  will  lift  fluid  away  from  the  wall  are  to  be  associated 
with  the  low-speed  streaks  and  where  the  vortices  move  fluid  towards  the  wall,  this 
region  will  be  associated  with  the  high-speed  streaks.  It  is  noted  that  this  idealized 
model  only  corresponds  to  ensemble  conditional  averaging  of  YSL  structures,  and 
doe^  not  take  into  account  that  the  vortices  lie  along  the  wall  for  short  distances. 
Whatever  the  real  picture. vortical  structure  is  dominant  in  the  viscous  sublayer. 

A  two  eddy  component  model  is  chosen  for  initial  development.  The  two 
components  model  the  small  scale  eddies  1SSE)  together  with  organized  large  scale 
eddies  il.SE).  These  organized  large  scale  eddies  exist  on  top  of  the  YSL  and  have 
length  scales  of  order  b.  Turbulence,  of  course,  exhibits  a  wide  range  of  scales  in 
reality.  This  may  be  seen  from  figure  10  showing  spectral  density  data  of  Fulachier 
( 1972)  for  the  t hree  velocity  components  and  the  temperature.  It  is  to  be  noted  that 
the  temperature  exhibits  scales  similar  to  the  streamwise  velocity  component.  The 
attempt  here  is  try  to  model  the  main  characteristics  of  the  viscous  sublayer  with 
only  two  eddy  scales;  one  scale  X:  corresponding  to  the  spanwise  spacing  between 
streaks.  and  another  scale  (  much  larger  than  A.  corresponding  to  boundary  layer 


thickness.  The  subscript  e  refers  to  outer  edge  values  ( y~  =  40).  The  additional 
subscript  1  refers  to  the  SSE  component  and  the  subcript  2  refers  to  the  LSE 
component.  The  general  form  of  the  boundary  conditions  applied  at  y~  -  40  is 


Component]  -  SSE  Component  2  -  LSE 
scale  A,  scale  6 


H  C  =  u,i(j.ye.z.f)  - 

«e2(0 

rf  -  rf  \  (x.yt.  z.t)  - 

r«2(0 

u~  -  ire  \(i.yc.  z.t) 

u'c2(0 

(3.1) 


0  -  6  -  6~  -  yt.  z.t)  -  0t2[t) 

The  SSE  which  are  associated  with  the  streaky  vortical  structure  are  respon¬ 
sible  for  the  production  of  turbulence  and  Reynolds  stress.  Visual  flow  observations 
indicate  that  the  SSE  can  have  characteristic  length  scales  as  small  as  10-30  wall 
units.  The  LSE  on  the  other  hand  do  not  produce  Reynolds  stress.  They  exist  in 
the  region  external  to  the  YSL.and  must  be  modeled  because  of  their  interaction 
with  the  SSE  and  their  contribution  to  pressure  fluctuations.  The  LSE  are  modeled 
an  unsteadv  stream  oscillating  at  the  outer  edge  of  the  YSL. 

This  two  component  model  is  compatible  with  the  concept  of  “active"  and 
“inactive"  components  of  turbulent  motion  in  the  log  region  as  characterized  by 
Townsend  (1901)  and  Bradshaw  (1967).  The  SSE  component  is  active  producing  the 
Reynolds  stress,  rotational,  and  dependent  on  wall  variables.  The  LSE  component 
is  inactive,  producing  pressure  fluctuations,  irrotational.  and  dependent  on  outer 
variables.  There  is  much  experimental  information  that  turbulent  flow  in  the  log 
region  and  hence  at  the  outer  edge  of  the  viscous  sublayer  comprises  these  two 
distinct  types  of  motion.  However.  Coles  (1978)  concluded  that  a  three  component 
model  is  needed  to  describe  the  turbulent  dynamics  in  the  viscous  sublayer.  The 
third  component  h  needed  a1-  a  t ransit iona)  component  between  the  “active  SSE 


and  the  "inactive"  LSE.  so  the  use  of  only  two  components  may  not  be  sufficient  1\ 
realistic  for  some  turbulence  parameters.  Other  aspects  and  complexities  of  this 
type  could  be  added  without  difficulty  in  principle. 

In  ari  effort  to  incorporate  coherent  structure  observations  in  a  relatively 
simple  two  component  model,  the  boundary  conditions  at  the  outer  edge  of  the 
viscous  sublayer  (y~  ~  40  ).  as  expressed  in  conventional  nondimensional  wall 
variables,  are  structured  as  terms  harmonic  in  time  and  span  as  follows: 

Component  1  -  SSE  Component  2  -  LSE 
scale  A  ,  scale  t 

'p  -  2o  ]  sin  .Vj  Tsin  c  -  ^  2(o-  -  o*)  s\n{.\\2T  *  Ou2) 

(3.2) 

v,  -  2J\  sin(  A  1 7  -  ox  i  j  sin  ,* 

wt‘  --  2oj  sin(.V)7'  -  o„j)cosc  -  ^  2 (a-  -  o* )  sin(.Vu27  -  ou.:>) 

0,'  2i]  sin(.V)7’  -  Cvjlsinc  -  ^  2(v*  -  if )  si  n  (  ACT  -  oP2) 

w  i  t  h 


A-  (Spacing 

between  streaks) 

(3.3  a) 

; 

2~  ~ 

A  - 

(3-3 b) 

a  ip  .1  i  . 

a  -  v[  v  -  e[ 

(3.3 c) 

\  <  >  ]  <  i  r  i  i  j  r> 

\ c)  -  <•$  v^\  vi  *  ’>'2 

(3.3  d) 
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O  1 

cos(o,  ]  ) 
a 

(3.3c) 

(/c  .  t,  1 

V 

1  COS(o,  1  -  Oh]  ) 

(3.3/) 

:*  cos (o-p.  Ou2  -  a  )  v i  ! 

\  h  — 

ot 

(3. ;•;<?) 

Experimental  measurements  give  various  information  about  the  boundan 
rnndit  ion  parameters  in  t  he  cornputat  ional  model.  The  low  speed  streaks  seen  in  t  he 
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viscous  sublayer  have  a  spacing  of  about  X~  =r  100.  This  gives  the  spanwise  length 
of  the  computational  domain.  The  vortical  structure  seen  by  experimentalists  is  not 
really  periodic,  however  the  vortical  pairs  have  a  measurable  mean  frequencv  and 
occur  repetitively  enough  in  the  viscous  sublayer  that  the  term  “quasi-periodic" 
applies  (  Blackwelder  and  Eckelmann  1979  ).  In  the  present  model,  the  streaky 
structure  is  assumed  to  be  periodic,  hence  periodic  boundary  conditions  are  applied 
spanwise.  The  streaks  which  alternate  with  their  high  and  low  speed  parcels  of  fluids 
provide  a  guide  for  modeling  the  spanwise  variation  of  the  streamwise  velocity  it 

The  coherent  streamwise  velocity  structures  are  highly  elongated  in  the 
streamwise  directions  with  X~  =  1000.  This  leads  to  a  mathematical  simplification 
of  neglecting  derivatives  in  the  streamwise  direction  (x)  compared  to  the  derivatives 
in  the  spanwise  (z)  and  normal  (y)  directions.  This  basic  approximation  which  in¬ 
volves  the  computation  of  three  velocity  components  in  two  space  directions  has 
been  either  termed  “2.5D  flow",  “slender  turbulence  theory",  or  “slender  eddy  the¬ 
ory".  As  noted  before,  the  dominant  eddy  stuctures  in  the  viscous  sublayer  are 
streamwise  contra-rotating  vortices.  This  information  provides  a  guide  for  struc¬ 
turing  the  spanwise  variations  of  r',  and  u  ",.  To  account  for  the  contra-rotating 
\ortical  motion,  with  •;  2r_  A ,.  the  t  ',  term  is  given  a  -  sine  dependence  and 

the  ir,’,  term  is  given  a  cos.-  dependence.  The  u',  term  is  given  a  sine  dependence 
in  view  of  the  observations  of  conditional  samples  of  the  bursting  process. The 
term  i-  1M)  degrees  out  of  phase  with  (Wallace  et .  al.  1977).  For  a  heated  wall. 
(?/,  is  in  phase  with  r,',.  A  conditional  sample  record  of  velocity  and  lemperature 
in  a  bursting  process  is  shown  in  figure  11  to  illustrate  this  point  (Antonia  et.  al. 
1982).  This  follows  the  conventional  picture.  In  the  ejection  process  for  a  heated 
wall,  a  temperature  parcel  transported  up  (?•  •  0)  will  be  at  a  higher  temperature 
than  its  time-mean  level.  Therefore  i.O.  and  v6  will  all  be  positive. 

An  essential  constraint  must  be  placed  on  in  order  to  conserve  mass. 
Since  periodic  boundarx  conditions  arc  used  in  the  spanwise  direction,  there  is  no 


1  1 


net  mass  flux  in  the  spanwise  direction.  Also  at  the  wall  no  mass  can  flow  through. 
At  the  outer  boundary,  therefore,  it  is  required  for  continuity  that 


(3.4) 


Because  tyo  is  a  function  of  time  alone,  it  must  be  set  to  zero  for  global  mass 
continuity  at  every  instant  of  time.  Additional  components  for  ve  obviously  could 
be  added  provided  such  components  have  a  spatial  dependence  which  satisfies  the 
global  mass  continuity  requirement. 

In  modeling  the  time  variations  of  the  velocity  and  temperature  boundary 
conditions,  sinusoidal  functions  are  used  with  allowance  for  phase  and  frequency 
differences  between  the  fluctuating  components.  The  frequency  is  set  by  the  mea¬ 
sured  frequenc\  of  both  ejection  sweep  and  sweep,  ejection  events.  Both  events 
occur  equalh  as  often  in  the  limit  of  low  threshold  and  long  integration  time  as 
shown  by  Johansson  and  Alfredson  (1982).  The  structure  of  the  present  boundary 
conditions  comprises  equal  proportions  of  both.  A  representative  burst  frequency 
has  been  found  from  the  experimental  measurements  of  Blackwelder  and  Haritoni- 
dis  1 1982),  Their  mean  frequency  of  the  ejection  sweep  event  is  .0035  averaged 
o\.  r  space  and  time  in  wall  units.  The  frequency  for  both  ejection  sweep  and 
sweep  ejection  events  will  be  twice  that  using  the  observation  of  Johansson  and 
Alfredson.  For  radian  measure,  a  factor  of  2 7t  needs  to  be  added.  Using  this  infor¬ 
mation.  the  bursting  frequency  .Yj  has  a  value  of  .044.  The  mean  period  of  the  LSE 
is  taken  as  T i > y  =  bf>  f'-v  (Laufer  and  Narayanan  1971).  This  sets  the  frequency 
Y,  -  Tise  2tt.  The  Reynolds  number  has  been  chosen  so  that  Re?  is  near  103. 
Measurements  of  T p  >  in  this  Re  range  have  been  made,  and  this  avoids  the  prob¬ 
lem  of  whether  T p  >  scales  with  wall  variables  or  free  stream  variables.  For  this 
Rt~.  it  turns  out  that  A',, :  A  ,.  The  values  of  the  intensities  of  the  velocities  and 
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the  temperature  are  also  used  to  set  some  of  the  parameters.  The  values  at  y~  -  40 
are  taken  as  u~'  =  a  =  2.0, 1'T*  =  3  =  1.0. u^1  =  a  =  1.3.  and  B~'  =  v  =  1.6. 

Using  appropriate  correlation  data  at  y *  —  40  most  of  the  coefficients  and 
phase  angles  can  be  determined.  Conditional  sample  data  of  the  burst  events  show 
a  180  degree  phase  difference  in  ucl  and  rel  which  sets  0vl  —  0.  Experimental  data 
for  a  heated  plate  also  show  that  uc;  and  #ei  are  also  180  out  of  phase.  This  sets 
or i  -  tt.  Correlation  coefficient  measurements  of  Ruv  -  -  .45  sets  qj  =  ,45a.  and 
of  R  -  -  .45  sets  L'  i  --  .45v.  The  Ru?  correlation  gives  a  relationship  between  ou 2 
and  o~-. 


Or 2  =  ou ■>  -  arccos( 


ctvRu *  - 


ft2t'2 


(3.5) 


The  correlation  Rut  has  a  value  of  -0.8  at  y~  =  40.  Also. 0^1  has  been 
fixed  in  order  to  be  consistent  with  experimental  data  showing  that  —  ■- -  %  0. 


dy 


From  the  boundart  conditions  of  the  mode). 


d[v;')- 

dy 


4rOij) 

c  os  ( o  u  1 

A  - 


o,  )  ) 


(3.6) 


Since  o  1  0  previously,  then  ou  j  -  *  2. 

This  still  leaves  some  remaining  parameters  which  have  not  been  fixed  by 
experimental  data.  These  are  a  j .  A  u  2.  o,.  2.  ©u2 .  and  A'r2.  Three  of  these  are  set  by 
assumption  and  and  two  by  computational  experimentation.  Thus,  cj  is  set  equal 
to  Jj.  since  t fie  vortical  structure  should  have  comparable  magnitudes  for  vjt  and 
u  ~  .  Using  the  assumption  of  staggered  LSE.A'u2  has  been  set  to  .Vu2  2.  In  view 
of  the  similar  spectral  energy  distributions  of  6  and  u.  AU2  is  set  equal  to  .Vu2.  The 
two  phase  angles  set  by  computational  trial  are  ou2  -  1.05  and  ou  2  -  3.14  radians. 
These  parameters  have  slight  effects  on  the  mean  velocity  profile  u(y).  u'.and  uv 
and  were  selected  to  best  fit  the  law  of  the  wall  for  u(y)  and  0(y). 
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The  use  of  sinusoidal  functions  for  the  time  and  space  variations  as  bound¬ 
ary  conditions  for  the  viscous  sublayer  model  would  fall  into  the  category  of  “psuedo- 
turbulence" .  This  type  of  analysis  has  been  useful  in  turbulence  analysis  especially 
in  allowing  analytic  solutions  to  be  obtained  for  certain  idealized  flows.  These  are 
not  real  turbulent  solutions,  but  they  provide  some  insight  about  turbulent  flow- 
(Hinze  1975}.  Most  of  these  analytic  attempts  reduce  the  Navier-Stokes  equations 
to  some  linearized  form.  The  present  investigation  with  boundary  conditions  which 
are  sinusoidal  functions  differs  slightly  since  the  full  Navier-Stokes  equations  are 
used  for  the  computation. 

A  not  her  point  about  t  he  use  of  sinusoidal  functions  for  boundary  condit  ions 
concerns  the  bursting  event.  I  se  of  sinusoidal  functions  in  the  two-component 
model  does  not  allow  for  any  negative  interactions  (u  >  0  and  v  >  0  or  a  <  0  and 
r  -  0).  During  the  burst,  these  have  been  measured  near  the  outer  egde  of  the  YSL 
to  contribute  negatively  about  .°>0(7  to  the  Reynolds  stress  (Wallace  1972).  The 
ejection  events  (t;  •  0  and  r  0 )  and  sweep  events  ( u  >  0  and  r  -  01  near  the 

outer  edge  of  the  VSL  contribute  about  70(7  and  GOT  respectively  to  the  Reynolds 
stro«s.  thi-  model,  there  are  no  negative  interactions  and  the  ejection  and  sweep 
events  contribute  equally  to  the  Reynolds  stress  with  50(7  each.  Also  the  burst 
process  jv  a  highly  intermittent  process  in  time  and  space.  The  bursting  process 
makes  up  only  about  20  25c7  of  the  mean  period  between  bursts.  W'hen  the  burst 

i-  not  occurring,  a  relativley  quiescent  signal  is  observed.  This  intermittency  is 
tiot  being  modelled  with  the  current  sinusoidal  functions.  A  Fourier  series  type 
boundary  condition  is  one  way  the  boundary  conditions  can  be  modified  to  model 
t  he  intermittency  better. 

With  the  abewe  considerations,  the  velocity  and  temperature  boundary  con¬ 
ditions  for  the  two-component  coherent -st ructure  model  of  viscous  sublayer  turbu¬ 
lence  become: 
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Component  1  -  SSE  Component  2  -  LSE 


scale  A,  scale  6 

u  ~  ~  2o  !  sin  .V]  T  sin  <■  -  ^  2(a:  -  q*)  sin(  A'u2T  -  ou2) 

(3.7) 

i\  =  -  23\ sin ( A  1T)  sin  c 

u  ~  -  2cr i  cos(A‘1T)  cos  c  —  ^  2(<72  -  o sin(A’u27'  ~  0^2) 

0;  -  -2fi  sinlA'iE)  sin  c  -  y  2(b’2  -  vf)  sin^V^r  -  092) 

with 

A'i  -  A'u2  -  A ~o  =  -044  (3.8a) 


ou2  =  1.05  ou2  =  3.14  or2  =  3.47  (3.86) 

o  =  2  a,  =  .9  J  =  J,  =  1.0  o=1.3  Oj  =  1.0  v  -  1.6  v,  =  .72  (3.8c) 

These  boundary  conditions  can  now  be  applied  to  a  time-dependent  Navier-Stokes 
code  and  turbulent  quantities  can  be  computed.  Figure  12  shows  the  computational 
domain  and  its  relation  to  the  boundary  layer.  When  these  boundary  conditions 
are  applied  to  the  Navier-Stokes  equations,  all  normal  (y)  derivative  quantities  are 
not  fixed.  These  normal  derivative  quantities  are  computed  as  part  of  the  numerical 
solution  of  the  Navire-St  okes  equations.  Therefore,  for  each  time  step,  the  boundary 
condition^  with  the  computed  normal  derivative  terms  will  be  consistent  with  the 
Navier-Stokes  equations  at  the  boundaries.  The  LSE  component  of  the  boundary 
conditions  are  being  modelled  as  an  unsteady  oscillating  flow.  To  be  consistent  with 
oscillating  2-D  shear  flow,  a  pressure  body  force  term  is  added.  This  term  will  be 
described  in  detail  later. 
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4.  Development  of  Computational  Code 

The  numerical  scheme  selected  for  the  viscous  sublayer  model  is  an  adap¬ 
tation  and  modification  of  the  Pulliam  and  Sieger  (1980)  code.  This  code  is  based 

on  a  three  dimensional  implicit  factorized  algorithm  for  the  unsteady  compressible 
Navier-Stokes  equations  with  the  thin-layer  approximation.  The  code  was  written 
for  arbitrary  geometries  and  even  though  the  computational  domain  for  the  viscous 
sublayer  model  is  cartesian,  the  use  of  metrics  makes  clustering  conveniently  easy. 
The  clustering  will  be  discussed  later. 

The  Pnlliam-Steger  code  had  to  be  modified  to  adapt  it  to  the  viscous  sub¬ 
layer  model  in  seven  ways  as  follows:  (1)  the  grid  has  been  stretched  to  cluster 

points  near  the  wall  and  the  outer  edge:  (2)  the  flow  variables  have  been  made 

nondirnensional  in  the  conventional  wall  variables  used  for  turbulent  flows:  (3)  var¬ 
ious  necessary  viscous  terms  have  been  added  to  bring  the  Pulliam-Steger  “thin 
layer"  code  in  the  full  Navier-Stokes  form:  (4)  the  boundary  conditions  have  been 
restructured  to  be  appropriate  fot  computations  of  viscous  sublayer  turbulence:  (5) 
the  algorithm  for  time  derivatives  has  been  revised  for  increased  accuracy  and  sta¬ 
bility:  ( f v j  the  dissipation  function  has  been  removed:  (7)  body  force  term  associated 
with  the  YSL  model  has  been  added:  (8)  f  differencing  has  been  dropped  using  the 
slender-edd\  approximation.  After  these  modifications  are  described,  numerical 
accuracy  and  stability  will  be  discussed. 

'4.1  Basic  Code  and  Finite  difference  algorithm 

\ariables  have  been  modified  slightly  to  conform  with  turbulence  conven¬ 
tion.  The  definition  of  the  physical  coordinates  and  velocities  in  the  Pulliam-Steger 
cod<  is  different  from  turbulent  boundary-layer  convention.  In  boundary  layer  con¬ 
vention.  the  streamwise  direction  and  velocity  are  x  and  u.  the  normal  direction  and 
velocity  are  y  and  v.  and  the  spanwise  direction  and  velocity  are  z  and  w. respect  ively. 
Pulliam  and  Steger  employ  a  different  coordinate  system  in  which  their  normal  di¬ 
rection  and  velociu  are  z  and  w,  and  their  spanwise  direction  and  velocity  are  > 
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and  v.  This  is  NASA  wing  convention.  This  convention  is  also  used  in  meteorology. 
In  order  to  minimize  code  changes,  their  coordinate  system  is  used  in  the  numeri¬ 


cal  computations  and  then  for  the  output,  the  boundary  layer  convention  is  used. 
Overbars  will  be  used  to  denote  the  Pulliam  and  Steger  variables.  Therefore,  the 
coordinates  are  x.y  =  z.:  =  y:  and  the  respective  velocities  are  u.r  =  u  .ir  -  r. 

The  grid  is  basic  with  clustering  only  in  one  coordinate.  Simple  cartesian 
coordinates  are  used  for  the  streamwise  and  spanwise  directions.  A  stretched  mesh 
is  used  in  the  normal  direction.  Hyperbolic  tangent  clustering  is  used  so  that 
points  can  be  concentrated  both  near  the  wall  and  the  outer  edge.  This  allows 
for  greater  resolution  near  the  boundaries.  The  metrics  for  this  transformation  are 
].yc  -  l  iy'.y,  =  1  Ac*  where  A y~  and  A z~  are  the  mesh  spacings 
in  wall  units  in  the  spanwise  and  normal  directions  respectively.  The  stretching 
function  used  is  the  following 


tanh{ ! 
a 


2  [L 

(/max 


1 1  tanh  1  ( a ) } ) 


(4.1a) 


c:  (4.16) 

nun 

with  a  0.9s  and  -  number  of  grid  points  in  the  normal  direction. 

With  the  grid  and  Pulliam-Steger  notation  defined,  the  nondimensional 
strong  conservation  form  the  the  Navier-Stokes  equations  with  all  the  approxima¬ 
tions  of  the  model  included  in  wall  variables  is  now  presented. 
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The  interna)  energy  is 
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0.5(u"2  -  v~  2  -  ir~2) 


and  pressure 


p  -  -  1)  e~  -  O.op~  (u~2  -  i  ~ 2  -  U'~‘)  (4.9) 

Temperature  is  related  to  pressure  and  density  using  the  equation  of  state  for  a 
perfect  gas 


9"  =  p ' 

A  lr  --Ur  \  ^R6  f 

where  the  friction  temperature  6  <  is  defined  by 


(-4.10a) 

(4.10&) 


Conventional  wall  variables  are 


pc.u- 


(4.10c) 
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V  'u  Pv  P  = 


6~  =  7T  p~  = 


U 

V  _  U' 

u~  =  — - 

Hr  UT 


T  =  r  = 


.X  =  n 


(4.11) 


In  this  code  0  represents  the  absolute  temperature  in  wall  variables.  When 

plotting  temperature,  a  modification  is  made  so  that  0~  -  0.  For  a  heated  wall. 

6U  ~  6 

this  becomes:  6~  —  ( — - - ). 

V  ( 

J 

The  finite-difference  scheme  used  is  the  Beam- Warming  (1978)  approximate 
factorization  algorithm  in  delta  form. 


(I  -  hSriBn  -  h6riJ-lBlJ)(l  -  hS-Cn  -  h6:J~xCvJ){qn *'  -  qr 
=  ~htv(Fn  -  F,n)  -  hb-(Gn  -  6’" )  -  hFr  -  .56A<?n 


(4.12a) 


a ,  2A' 

3 


(4.126) 


fi.B,  .C.C\  are  the  Jacobian  matrices  obtained  from  time  linearization  of  F  .FV.G,G  v 
respectively.  These  are  written  out  explicitly  in  Appendix  1.  The  fourth  order 
dissipation  operator  typically  used  has  been  removed.  Various  modifications  to 
the  code  were  necessary  to  make  it  compatible  with  the  present  computational 
model.  Differencing  in  the  J  direction  has  been  dropped  because  of  the  slender- 
eddy  approximation.  The  addition  of  the  complete  set  of  viscous  terms  is  needed 
since  spanwise  motions  are  to  be  modeled  by  the  Navier-Stokes  equations.  Spanwise 
shear  terms  were  not  included  in  the  Pulliam-Steger  code  because  of  their  thin- 
layer  approximation.  The  time-differencing  was  adapted  to  a  three-point  backward 
scheme  at  the  suggestion  of  Dr.  T.  Pulliam.  The  code  modification  was  minimal  for 
this  time  algorithm  change.  This  modification  was  made  for  reasons  of  stability  as 
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shown  by  Beam  and  Warming.  The  nondimensionalization  was  also  changed  from 
free-stream  reference  used  by  Pulliam  and  Steger  to  the  conventional  turbulent  wall 
variables  as  defined  previously. 

The  term  not  common  to  the  Navier-Stokes  equations  is  the  body  force 


term  F.  The  addition  of  the  right  hand  side  term  of  Fp  allows  for  a  forcing 
func  tion  compatible  with  oscillating  shear  flow.  The  function  Fr  plays  the  role  of 

the  pressure  fluctuations  caused  by  the  LSE  interaction  with  the  viscous  sublayer 

,  dp(  ,  dpt  dwe 2  ,  •  j 

b\  setting  -  ~  p- -  and  — —  =  - p  — — •  W  ith  these  as  the  imposed  pressure 

dx  dt  az  at 

gradients,  the  model  is  compatible  with  the  “active  inactive"  concept  mentioned 


pre\  lou'h 


In  the  VSL  model,  there  is  no  reason  to  assume  that  the  external 


•'lire  gradients  are  correlated  one  to  one  with  the  LSE  components.  It  has  been 


found  that  when  the  body  force  terms  were  not  used,  the  phasing  of  the  ur  and 


term'  near  wall.y*  -  2.  was  opposite  in  sign  to  experimental  measurement. 

Therefore,  the  body  force  terms  were  necessary  to  get  the  near  wall  phasing  of 


these  two  components  consistent  with  experiment.  Also  when  a  mean  streamwise 

procure  gradient  i'  applied,  a  constant  corresponding  to  the  pressure  gradient  is 
dp  ' 

added  to  the  -  term  in  F... 
dx  * 

12  C  ode  Verification  and  Accuracy  Test 


Since  the  Pulliam-Steger  code  is  a  compressible  code,  there  is  some  concern 
about  accuracy  and  stability  since  the  computation  is  done  at  low  subsonic  Mach 
numbers.  Ising  a  compressible  code  for  this  computation  is  not  necessary,  however, 
computations  of  this  type  increase  the  genera!  understanding  of  how  these  codes 
perform  in  the  limit  of  incompressible  flow.  An  empirically  set  lower  bound  to 
guarantee  numerical  accuracy  and  stability  is  about  M  >  0.2.  The  code  for  this 
research  is  being  run  in  a  low  Mach  number  range  (about  .1).  Therefore  since  the 
code  for  this  computation  is  being  used  in  a  very  extreme  case,  a  test  flow  must 
be  computed  to  show  that  the  numerical  scheme  is  accurate  and  stable.  Oscillating 
shear  flow  with  isothermal  and  adiabatic  temperature  boundary  <~onditions  was  used 


as  the  test  flow  for  assessing  the  accuracy  of  the  time  derivative  algorithm  and  the 
stability  of  the  present  code  at  low  Mach  numbers.  The  boundary  conditions  for 
oscillating  shear  flow  are 


dp  du  e 

i  -  u  -  ut  -  A  stn  nt  —  =  -p  — 

ox  at 


u  e  —  B  sin {mt  —  o) 

6f  =  C 


=  -P- 


(4.13) 


with  u  -  Dy.  and  A.B.C.and  D  are  arbitrary  constants.  The  analytic  solution 
for  the  velocity  fields  is  given  b\  Chapman  and  Kuhn  (1981).  The  temperature 
solution  has  been  worked  out  and  is  presented  in  Appendix  II.  Test  runs  of  the  code 
for  oscillating  shear  flow  gave  very  accurate  results  as  can  be  seen  in  figures  13  and 

Q 

14.  These  test  runs  are  for  a  heating  ratio  =  1 .0(adiabatic)  and  l.l(heated). 
Mach  number  \1  =  .02.  and  Pr  =  1.0.  This  test  show-s  that  time  accuracy  is  good 
and  the  compressible  code  is  stable  for  low  Mach  numbers.  These  results  show  that 
the  errors  due  to  the  cross  terms  from  the  approximate  factorization  do  not  affect 
the  time  accuracy.  It  was  found  for  stability,  however. that  the  parameter  A/r  from 
the  equation  of  state  should  be  less  than  0.4. 

14.3  Decoupling  and  Stability 

In  running  the  algorithm  with  the  desired  boundarv  conditions,  several 
numerical  problems  were  encountered  in  the  course  of  research.  A  brief  description 
of  these,  and  the  methods  devised  to  solve  or  minimize  them  are  described  in  the 
following  paragraphs. 

Global  mass  continuity  needs  to  be  satisfied  or  else  there  is  a  density  and 
pressure  build-up.  In  order  to  guarantee  global  mass  continuity,  it  was  found  that 
the  specification  of  mass  flux  variables  at  the  boundaries  was  necessary.  When 
just  the  velocities  were  specified,  a  large  build-up  of  density  and  pressure  was  ob¬ 
served.  Consequently,  the  velocity  boundary  conditions  were  changed  to  boundary 


conditions  on  pu~  —  p(ue,pv~  =  pfie.  and  pu\  =  pewe.  This  corrected  the  prob¬ 
lem  of  the  severe  density  and  pressure  build-up.  however,  other  problems  with  the 
boundary  conditions  were  encountered. 

Imposition  of  boundary  conditions  on  a  compressible  algorithm  still  does 
not  have  precise  rules  to  follow  and  thus  numerical  difficulties  can  arise  if  incorrect 
boundary  conditions  are  imposed.  Initially,  explicit  boundary  conditions  were  im¬ 
plemented  in  the  code.  At  the  outer  edge  of  the  viscous  sublayer,  four  boundary 
conditions  have  been  constructed  (equations  3.7).  At  the  wall,  the  no  slip  velocity 
boundar>  conditions  and  an  isothermal  temperature  boundary  condition  also  give 
four  boundan  conditions.  Five  boundan  condition?  are  required,  however,  for  a 
three  dimensional  compressible  code.  Conventional  numerical  schemes  use  the  en- 
erg\  equation  or  the  normal  momentum  equation  to  back  out  pressure  as  the  fifth 
boundan  condition  at  walls.  For  the  spanwise  direction,  periodic  conditions  were 
irnpo-ed  on  all  five  variables  a^  stated  before.  This  was  implemented  by  using  a 
block  periodic  sober  Initially  a  normal  momentum  equation  is  used  as  the  fifth 
boundan  condition  to  back  out  pressure  both  at  the  wall  and  the  outer  edge. 

With  these-  <  omen:  ional  explicit  boundan  conditions,  the  code  was  either 
unstable  or  tn«  solution  exhibited  decoupling.  Decoupling  is  when  the  solution 
on  adjacent  grid  points  is  discontinuous,  in  a  saw-tooth  fashion,  but  solutions  on 
even  other  grid  point  are  continuous.  Therefore,  there  exists  a  smooth  solution 
for  the  cut.  points  and  another  for  the  odd  points  when  decoupling  occurs.  This 
is  a  t\pical  problem  when  central  differencing  schemes  are  used.  Decoupling  of 
a  solution  is  often  a  signal  that  there  is  a  problem  with  the  boundary  condition 
specification.  A  signal  that  a  problem  exists  at  the  outer  boundary  can  also  be  seen 
from  the  relatively  high  mass  continuity  residual  which  exists  there.  At  the  outer 
edge,  the  mass  continuity  residual  is  about  .01 .  whereas  interior  continuity  residuals 
are  in  the  range  of  10  f  The  artificial  boundary  conditions  cause  a  difficulty  for 
the  continuity  equation  to  be  satisfied  numerically  at  the  outer  edge. 
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Problems  with  decoupling  were  most  noticeable  when  the  oscillating  normal 
velocity  component  i\  was  nonzero.  Some  of  the  problems  are  assumed  to  be 
associated  with  the  pressure  pulses  generated  by  the  oscillating  vf  component.  The 
decoupling  mainly  was  observed  in  p.p'.v.v'.  and  e.  Decoupling  was  most  noticible 
at  the  walls.  Also,  decoupling  became  more  severe  as  Prandtl  number  was  increased. 

the  heating  ratio  --  decreased,  or  the  friction  velocity  increased. 

0 1 

A  common  method  for  solving  the  decoupling  problem  is  the  use  of  an  ar¬ 
tificial  dissipation  function  The  fourth  order  artificial  dissipation  function  appear¬ 
ing  in  the  Pulliam-Steger  code,  however,  has  been  turned  off  for  ail  computations 
presented  here.  The  artificial  dissipation  function  is  commonly  used  in  central- 
differencing  algorithms  like  the  Pulliam-Steger  code  to  prevent  decoupling  and  to 
damp  high  frequencies  due  to  nonlinearit ies  of  the  Navier-Stokes  equations.  The 
artificial  dissipation  function  allows  for  more  coupling  between  the  computational 
grid  points  which  tends  to  eliminate  decoupling  of  the  solution  which  can  occur 
with  central  difference  algorithms.  It  was  found  that  when  the  artificial  dissipation 
function  is  used  in  the  YS1.  computation. the  spatial  location  of  peaks  of  quantities 
such  a1-  u'  are  shifted  away  from  experimental  measurements.  The  experimental 
peak  of  u'  is  at  about  y'  12.  The  computations  to  be  presented  show  a  peak 


at  about  y 


13.  (See  figure  25)  When  dissipation  was  used  the  peak  location 


moved  out  beyond  a  ~  25.  if  the  dissipation  was  very  large,  then  the  peak 
was  wiped  out  completely.  Also  dissipation  was  tried  only  in  the  energy  equation 
to  try  to  damp  out  high  values  of  p'\  but  this  proved  to  be  unstable.  Localized 
dissipation  near  the  boundaries  was  also  investigated  since  the  production  of  the 
high  p'  occurs  at  the  boundaries.  This  was  also  unstable.  Therefore,  the  dissipation 
function  was  not  used  for  any  of  the  computations  of  sublayer  turbulence.  The 
decoupling  problem  was  not  solved  using  dissipation,  therefore,  the  boundary  con¬ 
ditions  werer  econstructed  to  eliminate  or  minimize  the  decoupling  problem.  The 


damping  of  high  frequencies  from  the  nonlinearities  is  hoped  to  be  provided  by  the 
natural  dissipation  of  the  viscous  terms. 

Alternative  boundary  conditions  were  tried  in  an  attempt  to  fix  the  decou¬ 
pling  problem.  A  characteristics  type  approach  was  used  to  more  accurately  follow 
the  direction  of  the  propagation  of  information. The  signal  direction  was  specified 
bv  the  sign  of  the  vertical  velocity  v~  .  The  boundary  conditions  for  v~  >  0  are  out¬ 
flow  for  u~  ,v~  .  u'~  .  p~ .  and  6~  specified.  Fort/  <.  0.  u~ .  v~  .  U'~  .  8~  were  specified 
with  outflow  for  p~  .  Outflow  is  simple  first  order  extrapolation.  This  set  of  bound¬ 
ary  conditions  is  more  physically  correct  for  a  compressible  code:  however,  these 
boundary  conditions  did  not  help  the  decoupling:  and  for  Prandtl  number  greater 
than  .72.  the  code  was  more  unstable.  Other  variations  of  this  type  of  boundary 
condition  were  attempted,  but.  none  of  them  helped  solve  the  decoupling  problem 
and  some  even  did  not  capture  the  vortical  dynamics  needed  for  the  present  model. 

The  characterist  ic  method  imposes  Neumann  type  (gradient )  boundary  con¬ 
dition;  when  outflow  boundary  conditions  are  needed.  These  boundary  conditions, 
however,  have  been  found  to  be  much  more  unstable  than  the  Dirich  let  boundary 
conditions  for  this  computation.  A  modified  boundary  condition  was  tried  to  see  if 
the  decoupling  problem  could  be  helped.  A  linear  mixed  type  boundary  condition 
of  the  form 


ait  -  b  -  =  a f  i  -  b}_ 

dy 


(A.l-1) 


was  implemented.  This  should  help  decoupling  since  the  gradient  term  will  couple 
the  point  next  to  the  outer  boundary  to  the  applied  boundary  condition  when  one¬ 
sided  differencing  is  used.  This  type  of  boundary  condition  proved  to  be  stable, 
however  it  did  not  help  any  of  the  numerical  problems  compared  to  the  Dircihlet 


boundary  conditions  and  therefore  was  not  used. 


In  view  of  these  results,  it  was  decided  to  stay  with  the  imposition  of  Dirich- 
let  boundary  conditions  for  the  velocity  and  the  temperature  at  the  outer  edge  of 


ihe  viscous  sublayer.  This  requires,  however,  that  the  fifth  boundary  condition  be 
reformulated  to  alleviate  the  numerical  problems  since  the  normal  momentum  equa¬ 
tion  initially  used  was  inadequate.  Note  that  the  viscous  terms  give  the  equation  set 
an  elliptic  character  which  may  allow  for  a  full  Dirichlet  boundary  condition  speci¬ 
fication  Th<  selection  of  the  fifth  boundary  condition  under  these  conditions  is  still 


unclear  and  must  be  chosen  by  computational  experimentation.  Many  boundary 
londitions  of  various  Dirichlet  Neumann  combinations  were  tried.  For  example, 
a  radialiw  type  (non-reflecting)  energy  boundary  condition  (Bayliss.  Turkel  1981) 
wa-  tried,  but  this  boundary  condition  made  the  code  unstable.  All  attempts  to 
rectify  the  decoupling  problem  in  this  manner  failed.  Many  combinations  resulted 
m  severe  instabihies.  It  was  found  that  with  explicit  boundary  conditions,  the 
decoupling  could  not  be  solved. 

Since  all  formulations  with  an  explicit  boundary  condition  failed,  an  implicit 
boundary  condition  formulation  was  implemented.  Implicit  boundary  conditions  are 
said  to  be  generally  more  stable  and  robust  by  Beam  (1981).  As  with  the  explicit 
formulation,  it  wa>  decided  to  go  w  ith  the  Dirichlet  specificat  ion  of  the  velocities  and 
temperature.  The  fifth  boundary  condition  which  greatly  minimized  the  decoupling 
problem  wa<  found  to  be  the  imposition  of  the  full  continuity  equation  at  the  wall 

and  the  outer  boundary.  This  was  implemented  by  expressing  the  time  derivative 

c)  ft  idpx'} 

using  three-point  backward  time  differencing,  the  normal  flux  term  — — -  using 
dt  dy 

( dp  u* ) 

first-order  one-sided  differencing  and  the  spanwise  flux  term  — - -  using  second¬ 

ed 

order  central-differencing.  All  profiles  except  for  e  were  then  found  to  be  smooth 
and  well-behaved,  although  e  showed  only  slight  decoupling  near  the  wall  of  1  part 
in  10  (>ee  figure  15).  Therefore  the  continuity  equation  as  the  fifth  boundary 
condition  is  used  instead  of  the  normal  momentum  equation  at  the  wall  and  the 
outer  boundary  to  minimize  the  numerical  problems. 
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With  the  decoupling  problem  essentially  solved,  some  problems  of  stability 
still  remained.  Limitations  on  the  range  of  values  of  the  Prandtl  number  were 
encountered  when  the  Pr  was  varied.  In  the  physical  world,  gases  do  not  have 
Pr  •  1.  Since  the  perfect  gas  law  is  being  used  for  the  equation  of  state,  this  limits 
in  reality  the  Pr  range  to  order  one.  However,  since  Pr  dependence  is  a  desired 
parameter  to  study,  the  Pr  was  varied  from  .72  to  6.0  to  explore  the  effect  of  Pr. 
Code  instability  was  encountered  for  Pr  higher  than  10. 0.  This  may  in  some  way  be 
associated  with  the  use  of  a  perfect  gas  law  together  with  Pr  values  corresponding 
to  liquids.  This  is  obviously  unphysical. 

Stability  bounds  also  exist  for  the  Courant  number.  There  always  seemed 
to  be  a  minimum  Courant  number  for  stability.  The  minimum  Courant  munber  is 
about  50. This  limits  the  minimum  time  step  allowed  for  a  given  grid.  Run  conditions 
have  been  chosen  to  be  within  all  stability  bounds.  The  Courant  number  did  not 
seem  to  have  a  maximum  bound:  however,  accuracy  was  very  bad  if  the  Courant 
number  was  too  large  (30000).  (Remember  that  these  Courant  numbers  are  in  wall 
variables.)  The  time  step  limitation  shows  that  the  stability  region  is  of  limited  size 
and  does  not  include  the  region  near  the  origin. 

By  computational  trial  runs,  the  basic  stability  ranges  of  this  code  have 
been  determined  as  discussed  previously  in  this  section.  In  summary,  a  basic  list  of 
these  conditions  is  given. 

(1)  M-  <  0.4 

(2)  Pr  <  6.0 

(3)  (CFL)mir,  >  50. 

(4)  (CFL)max<  30000. 

The  viscous  sublayer  model  implemented  in  the  Pulliam-Steger  Xavier- 
Stoke'  code  was  then  run  on  the  CRAY-XMP.  and  computational  results  obtained 
for  various  turbulence  quantities. 
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5.  Computational  Results  and  Discussion 

In  starting  the  YSL  computation,  initial  conditions  at  time  zero  are  re¬ 
quired.  Initial  profiles  for  u  and  0  were  imposed  which  are  linear  piecewise  con¬ 
tinuous  and  match  the  mean  slope  at  the  wall  and  the  mean  slope  of  the  law  of 
the  wall  at  y~  40.  The  density  is  set  so  that  the  initial  pressure  is  a  constant 
using  the  perfect  gas  law.  Initial  values  for  r  and  tc  are  set  to  zero.  The  boundary 
conditions  on  velocity  and  temperature  (equations  3.7)  are  slowly  turned  on  using 
an  exponential  growth  cycle  of  two  A'j  cycles.  Periodic  flow  conditions  developed  in 
about  six  .Y,  cycles.  The  code  seems  to  be  relatively  insensitive  to  how  the  model 
boundary  conditions  are  started  up.  An  important  requirement,  however,  is  to  set 
the  correct  mean  values  for  uf  and  0f  in  the  initial  conditions.  If  this  is  not  done, 
the  mean  profiles  will  not  be  correct  even  after  periodic  conditions  are  achieved. 
Computations  were  made  for  ten  .Y,  c\cles.  Averages  over  space  and  time  starting 
from  the  seventh  A’,  cycle  were  taken  to  provide  the  various  turbulence  statistics. 
Figures  16(a), (b).(c). and  (d)  show  contour  plots  of  the  vertical  velocity  component 

i  before  any  averaging  C  applied.  It  can  be  seen  that  the  contours  for  t  at  A’jT  =  - 

4 

and  ^  are  consistent  with  contra-rotating  vortical  motion.  This  illustrates  that  the 
imposed  outer  edge  boundan  condition^  are  producing  the  desired  vortical  motion 
in  the  computational  domain. 

The  typical  run  conditions  were 

(1)  360  time  steps  per  A  ,  cycle 

(2)  39  grid  points  in  the  normal  direction.  25  grid  points  in  the  spanwise 
direction 

(3)  ■-  300  degrees  Kelvin. 6f  —  290  degrees  Kelvin 

(4!  Re  -  23.300.  -  50  m  s. 
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Run  times  were  about  600  seconds  CPC  time  on  the  CRAY-XMP.  To  check 
flow  periodicity  and  code  stability,  some  longer  runs  were  made  with  up  to  4000 
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seconds  CPI"  time  with  the  grid  essentially  doubled  to  79  points  in  the  normal 
direction  and  39  points  in  the  spanwise  direction,  and  with  time  step  also  refined. 

Output  quantities  such  as  the  mean  velocities  u.v.  w.  mean  density  p.mean 
temperature  fl.mean  pressure  p.  velocity  intensities,  u'.  v'.w'.  density  intensity 
p'.  pressure  intensity  p' .  and  temperature  intensity  O'.  Reynolds  stresses.  uT.vw. 
vTr ,i'd ,u6.u 0. correlation  coefficients,  Rul  .Rve,RuP  and  additional  correlations  such 
as  flatness  and  skewness  are  computed  as  space  and  time  averages.  The  results  for 
these  hydrodynamic  and  thermodynamic  quantities  in  the  YSL  for  Pr  —  0.72  are 
presented  graphically,  (see  figures  19-30)  Before  discussing  these  results,  however, 
a  brief  outline  is  presented  of  two  irregularities  observed  in  the  computations:  near 
the  boundaries,  and  in  the  pressure  fluctuations. 

;. 5 . 1  Irregularities  Near  Outer  Edge  and  Wall  Boundary 

In  developing  the  computational  model  it  has  been  found  that  control  of 
the  gradients  of  quantities  near  the  outer  edge  is  very  difficult  with  the  boundary 
condition  parameters  available.  Variation  of  these  parameters  seem  to  effect  the 
values  most  significant  1\  in  the  region  y~  <  35.  From  y~  =  35  -  40.  there  seems 
to  be  a  "btokes-type  lacer”.  where  the  physical  quantities  are  adjusting  to  the 
artificially  forced  motion  of  the  applied  boundary  conditions.  This  can  been  seen  in 
figure  17  which  shows  the  mean  gradients  of  streamwise  velocity  and  temperature. 
Anomalous  variations  exist  just  inside  the  outer  boundary  ( y~  ~  40)  .but  relax 
to  realistic  variations  at  about  a  y~  of  35.  This  anomalous  region  extends  slightly 
further  towards  the  wall  as  Pr  increases  and  if  a  mean  favorable  streamwise  pressure 
gradient  is  imposed.  The  adverse  pressure  gradient  does  not  seem  to  extend  the 
anomalous  region  at  all.  In  these  situations  the  anomalous  region  will  extend  down 
to  about  y  '  s:  30.  Because  of  the  relaxation  behavior  of  gradient  type  quantities  in 
our  model,  and  because  Pr T  is  made  up  of  two  gradients,  values  of  the  PrT  will  be 
restricted  to  y~  ■  35.  since  the  values  for  35  ■  y "  •  40  are  not  physically  realistic. 
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Al!  profiles  are  smooth  except  for  some  of  the  correlations  like  uv  and  v6 
which  exhibit  a  wiggle  in  their  values  very  near  the  wall.  (y~  <  .4)  This  is  a 
numerical  problem  which  may  still  be  associated  with  a  boundary  condition  problem 
at  the  wall.  Also  slight  decoupling  exists  in  mean  total  energy  e  of  1  part  in  1CF'  as 
mentioned  before  (see  figure  15).  This  decouping  is  localized  near  the  wall  and  its 
region  of  influence  increases  as  Pr  increases.  The  final  boundary  conditions  for  the 
wall  are  the  best  that  could  be  achieved  for  the  present  model,  but  still  yield  some 
wall  boundary  irregularities. 

The  computed  turbulent  Prandtl  number  involves  two  of  the  correlations 
which  exhibits  these  "wiggles",  and  hence  is  affected  by  the  near-wall  numerical 
problems.  Figure  18  shows  the  near  wall  behavior  of  Prp  for  Pr  —  0.72  as  the  grid 
is  refined  by  a  factor  of  2.  The  grid  dependency  shows  up  for  the  region  y~  <  0.4. 
For  y~  0.4  the  results  are  the  same.  These  wiggles  can  be  eliminated  with  a 
coarser  grid:  however,  the  finer  grid  is  desired  for  increased  resolution  near  the  wall. 
Because  of  the  small  grid  dependence  of  the  computed  values  near  the  wall,  values 
plotted  in  the  figures  are  restricted  to  y~  :■  0.4. 

Major  effects  of  grid  changes  show  up  in  the  computation  of  the  pressure. 
Also  effected  in  a  lesser  degree  are  density  and  temperature.  The  velocity  compo¬ 
nents  are  unaffected.  When  the  grid  is  progressively  refined,  a  sharp  kink  in  pressure 
at  the  wall  develops.  This  is  another  signal  that  there  still  may  be  a  boundary  con¬ 
dition  problem  at  the  wall.  The  wiggles  of  ux  .vO  and  e  near  wall  are  also  thought  to 
be  parth  caused  by  the  artificially  high  pressure  intentsity  at  the  wall.  Runs  have 
been  made  where  the  pressure  field  fluctuations  have  been  manually  reduced  each 
time  step  by  equal  proportions  so  that  p'u  =;  5  (reasonable  experimental  value). 
When  this  has  been  done,  there  is  no  decoupling  in  c  and  no  wiggle  in  ur.  This 
shows  that  the  high  pressure  intensity  at  the  wall  may  contribute  to  the  numerical 
problem  at  the  wali.  This  method  is  not  used  since  the  pressure  field  damping  is 
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arbitrary,  and  has  no  physical  basis  behind  it.  The  results,  however,  illustrate  quite 
well  the  effect  of  pressure  on  the  near  wall  quantities. 

.5.2  Pressure  Fluctuations 

The  computed  intensities  of  pressure  fluctuations  in  wall  variables  varies 

from  12  to  22  for  Pr  from  .72  to  6.0  and  for  —  from  —.01  to  -.01  (See  fig- 

dr 

ures  29(c). 31(d). 33(d). 35(d). 38(b). -40(b)).  These  values  of  p'  are  much  higher  than 
experimental  measurements  of  p'uall  which  are  in  the  range  of  3  -  5.  However 
the  computed  values  are  close  to  the  high  end  of  a  theoretical  range  developed  by 
Kraichnan  which  is  p'u  12  (Hinze  1975).  The  high  values  of  p'  are  believed  to 
be  caused  by  pressure  waves  generated  by  boundary  reflections  of  energy.  At  the 
outer  boundary,  very  steep  gradients  of  the  velocity  field  also  contribute  to  the  high 
pressure  intensity  through  the  pressure  gradient  term  of  the  normal  momentum 
equation.  A  solution  to  this  would  appear  to  be  the  construction  of  a  non-reflecting 
boundary  condition  for  energy,  and  the  minimization  of  the  steep  gradients  of  the 
velocity  field  at  the  outer  boundary.  Tnis  was  attempted  as  discussed  in  the  pre¬ 
vious  chapter:  but.  only  the  decoupling  problem  in  the  vertical  velocity  field  was 
solved. 

It  is  known  that  high  values  of  p'  can  be  a  sign  that  there  is  a  numerical  or 
physical  problem.  High  values  of  p'  signaled  the  fact  that  mass  was  being  pumped 
into  the  computational  domain  when  the  velocities  uf.r t.u.  rather  than  the  mass 
fluxes  were  being  specified.  The  solution  to  this  problem  .as  noted  earlier,  was  the 
change  to  mass  flux  specification,  pv  f.  pvt.  pu\  at  the  outer  boundary.  It  has  been 
found  that  p'  in  the  present  computations  is  affected  by  grid  spacing,  time  step, 
boundary  conditions.  A/-,  and  Pr.  These  changes  in  p' .  however,  seem  to  have 
negligible  effect  on  the  other  thermodynamic  and  flow  variables.  When  parameters 
were  chosen  near  a  stability  bound  (for  example  when  a  time  step  was  getting  too 
small),  very  large  p'  and  9'  values  were  generated.  In  these  cases,  the  very  large 
values  of  p  affect  the  temperature  field. 


The  value  of  p'  can  be  checked  at  the  wall  to  make  sure  its  value  is  correct 
in  relation  to  the  other  thermodynamic  variables.  This  can  be  done  by  evaluating 
the  equation  of  state  for  a  perfect  gas  at  the  wall  with  an  isothermal  wall  boundary 
condition  on  temperature.  This  relation  for  pressure  intensity  is 


(p'lu 


(*)uV)u 

“-M* 


(5.1 ) 


The  following  table,  using  the  values  of  [0)w.  [p')u,.  and  \iT  from  the  computational 
runs,  '-hows  that  the  computed  values  of  p’  are  consistent  with  the  equation  of  state. 


Table  5.1 

Comparison  of  values  of  (p')u 


Pr 

u. 

(p')u 

(p')u 

Eqn. 

[?')u 

Comp. 

0.72 

510.3 

0.1573 

.00101 

16.02 

16.03 

1.50 

746.0 

0. 1  CM 

.00082 

13.01 

12.98 

3.0 

To  44 .0 

*  0  2160 

.00088 

13.95 

14.00 

6.0  ‘ 

15~10.li 

0.26  if- 

'  .00007 

15.38 

15.44 

It  was  also  found  that  high  value-  of  p'  ran  be  generated  if  A/r  <<  1.  This 

is  evident  by  refering  to  equation  (5.11  In  these  cases  the  high  values  of  p'  seem  to 

have  little  effect  on  the  other  variables.  Again,  parameters  for  the  present  data  runs 

are  such  that  p'  =  12  -  22.  it  is  noted  here  that  for  all  computational  runs  and  for 

all  boundary  conditions,  p  is  essential!)  constant  throughout  the  YSL.  This  agrees 

with  the  conventional  steadv  state  two-dimensional  boundary  layer  assumption  that 
dp 

~  =  ° 

dy 

£5.3  Velocities. Temperature,  ami  Turbulent  Prandtl  Number 

The  computational  model  has  been  constructed  so  that  the  four  quantities 

which  make  up  PrT.  namelv.  ui  .vB.  .  and  -  are  as  realistic  as  possible  for  the 

<hj  dy 

given  boundary  condition-.  Modei  parameters  were  set  up  as  described  earlier  in  an 


effort  to  have  the  values  for  u.  6.  and  ur  fit  as  well  as  possible  the  experimental  value? 

for  Pr  -  .72  (Air)  with  zero  pressure  gradient.  Most  of  the  available  experimental 

data  for  6.  6'.  Ru„.  and  Rvr.  are  for  Pr  =  .72.  This  Prandtl  number  allows  for  the 

most  complete  comparison  of  the  model's  results  involving  temperature  terms.  The 

computational  and  experimental  data  for  Pr  -  0.72  have  been  plotted  in  figure^ 

19-26.  Additional  data  for  the  YSL  computation  can  be  found  in  figures  27-30. 

du  36 

The  mean  gradients  — -  and  —  are  compared  with  experiment  indirectlv 

dydy 

b\  looking  at  the  profiles  for  u  and  0.  Both  the  computed  mean  strearmvise  velocity 

and  mean  temperature  profiles  follow  the  experimental  data  fairly  well,  (see  figures 

du  dO 

IS  and  23)  The  slopes  at  the  wall.  ( —  )uait  and  (  —  )uall-  arc  within  5cx  of  their 

dy  3y 

dxi 

theoretical  values  of  1.0  and  Pr  =  0.72.  respectivelv.  The  mean  gradients  —  and 

dy 

36 

—  are  realistic  up  to  about  5:  36. 

dy 

The  Reynolds  stress  ur  is  consistent  with  experimental  measurements.  (See 
figure  21)  Values  of  ur  seem  to  be  slightly  high  near  the  outer  boundary  when 
compared  to  the  experimental  data  shown.  Overall,  the  profile  for  the  Reynolds 
stress  compares  well  with  experiment. 

Experimental  data  for  vB  is  not  directly  available  near  a  wall,  but  is  in¬ 
directly  available  in  the  form  of  the  correlation  coefficient  (See  figure  26).  The 
correlation  coefficient  /?rs  =  has  been  experimentally  measured  to  be  about 


.45  at  a  y' 


40.  This  value  has  been  used  a?  a  boundarv  condition  for  the 


model.  The  experimental  values  of  Rvn  had  to  be  slightly  extrapolated.  The  clos¬ 
est  data  points  for  R,  *  were  from  a  y~  of  about  50  to  70.  The  two  experimental 
data  sets  used  to  get  extrapolated  values  of  Rvt.  decrease  as  the  wall  is  approached 
(50  <  y~  ■'  100).  If  this  trend  is  to  continue  in  the  viscous  sublayer.it  would  be 
inconsistent  with  the  computed  results.  The  behavior  of  Rvp  can  be  expected  to  be 
similar  to  that  of  R,.,  for  a  heated  wall  as  long  as  6  is  behaving  as  a  passive  scalar. 
(Refer  to  figures  22  and  26)  The  computed  curves  for  both  exhibit  similar  behavior. 
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|  Townsend  (1976)  states  that  in  the  YSL.  the  correlation  Rv?  has  an  experimental 

range  between  .50  -  .55.  The  computed  results  are  in  this  range.  Also,  the  curves 

!  for  uv  and  r 6.  like  Rux  and  /?,  *.  are  similar  as  might  be  expected  (See  figures  21  and 

28(d)).  As  an  additional  check,  the  correlation  coefficient  Ru *  can  be  seen  in  figure 

25  to  compare  well  with  experimental  measurements.  Therefore  all  four  quantities 

jj  which  make:  up  the  Ptj  seem  physically  reasonable  for  Pr  -  .72. 

i  Once  the  YSL  model  boundary  condition  parameters  0U2 ■  Oe 2. and  ou-2  "ere 

set  for  Pr  =  0.72  and  zero  pressure  gradient,  the  model  was  held  constant  and  the 

dp 

1  parameters  of  Pr  and  — -  were  varied.  Data  for  the  varving  Prandt)  number  and 

:  di 

1  pressure  gradient  runs  can  be  found  in  figures  31-40.  Only  computational  results  for 

quantities  involving  temperature  .pressure  and  density  are  shown  for  the  different 

Prandtl  number  runs.  This  is  because  the  hydrodynamic  results  are  essentially 

constant  as  Prandtl  number  varies.  The  Pr  range  is  from  .72  to  6.0,  and  the  mean 

dp 

streamwise  pressure  gradient  was  varied  from  -.01  to  -.01.  With  changing  — 

dx 

a  corresponding  change  in  the  bursting  frequency  .Yj  of  the  SSE  was  employed 
to  conform  to  the  experimental  data  by  Schraub  and  Kline  (1965).  Values  of  the 
Bursting  frequency  .Yj  with  pressure  gradient  are  tabulated  below. 

Table  5.2  ' 


Values  of  .Yj  with  Pressure  Gradient 


dp- 

dx  - 

-V, 

-0.01 

0.0358 

0.0 

0.0440 

-0.01 

0.0583 

The  turbulent  Prandtl  number  Ptj  for  each  of  the  different  runs  can  now 
be  compared.  The  computed  results  for  variations  in  molecular  Prandtl  number 
Pr  for  an  isothermal  wall  are  presented  in  figure  41.  The  results  have  similar 
characteristic^  with  several  existing  theories  for  Ptj  (refer  back  to  figures  3-6). 
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The  Prj  is  about  one  for  the  whole  Pr  range  examined,  except  near  the  wall.  In 
the  region  y~  <  3.  Prj  is  stongly  affected  by  Pr.  As  Pr  is  increased.  PrT  near  the 
wall  increases  almost  proportionally.  Of  the  theories  listed,  only  the  ones  by  Geshev 
(1978).  P.L.  Maksin  et.  al.  (1977).  and  L.C.  Thomas  (1981)  show  an  increase  of 
PrT  as  Pr  increases.  The  others  show  opposite  behavior  or  no  dependence  on 
Pr.  Geshev's  model  is  from  a  formal  Green's  function  solution  to  the  linearized 
momentum  and  energy  equations.  Thomas'  model  is  a  modified  surface  renewal 
model.  Maksin's  model  uses  the  transport  equation  for  vO.  All  of  these  models 
including  the  present  model  have  been  developed  from  very  different  viewpoints, 
however,  all  are  qualitatively  alike.  The  increase  of  Prj  for  the  YSL  model  is  much 
stronger  than  for  the  theories  of  Geshev.  Maksin  and  Thomas.  The  model  of  Geshev 
has  two  different  wall  temperature  boundary  conditions,  isothermal  and  adiabatic. 
Geshev's  result  for  the  adiabatic  wall  boundary  condition  indicates  that  PrT  goes 
to  zero  at  the  wall  for  all  Prandtl  numbers.  Geshev's  result  for  the  isothermal  wall 
boundary  condition  shows  Pry  increasing  at  the  wall  with  increasing  molecular 
Prandtl  number.  Finite  peak  values  of  Prj  are  attained  for  finite  Pr.  For  Pr  —>  oc. 
then  Prj  —  oc.  Geshev's  analysis  indicates  the  near  wall  behavior  of  Prj  to 
be  dependent  on  the  wall  boundary  condition.  The  wall  temperature  boundary 
condition  for  the  present  computation  is  isothermal.  The  computed  Prj  for  the 
present  model  is  in  qualitative  agreement  with  Geshev's  model  for  isothermal  wall 
boundary  conditions.  The  computational  models  of  Thomas  and  Maksin  are  also 
for  isothermal  wall  boundary  conditions.  Again,  these  computed  results  of  Prj  are 
in  qualitative  agreement  with  the  isothermal  result  of  Geshev.  For  other  turbulent 
Prandtl  number  models,  various  other  temperature  boundary  treatments  have  been 
employed.  For  example.  Cebeci  also  allows  for  isothermal  wall  boundary  conditions, 
but  in  an  indirect  way  by  using  damped  wall  functions  for  the  eddy  diffusivity 
This  is  analogous  to  the  damping  functions  developed  for  the  eddy  viscosity  <  . 
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The  high  values  for  Prj  near  the  wall  for  Pr  -  6.0  are  also  in  qualitative 
agreement  with  the  results  of  Simonek  (1983).  He  states  that  the  behavior  of  r 8 
near  the  wall  changes  from  a  third  power  behavior  with  y~  to  a  fourth  power  when 
Pr  is  greater  than  7.  This  also  will  cause  Prj  to  increase  as  the  wall  is  approached, 
since  this  will  give  PrT  a  —  behavior.  This  studv  does  not  have  anv  Pr  greater 

y 

than  7.  however  a  check  can  be  made  with  the  Pr  =  6.0  computation  since  it  it  is 
close  the  the  "transition"  Pr  of  7. 

The  power  law  log-log  plots  for  Pr  =  0.72  in  figures  43-45  show  behavior 
close  to.  but  not  precisely,  as  v1  —  y~.r'  ~  y~2,w'  ~  y~.6'  ~  y~.uv  ~  y~3 
and  v8  ~~  y-3.  In  figures  46-48.  the  correlation  v6  is  plotted  logarithmically  to 
show  the  near-wall  behavior  as  Pr  is  increased.  Figure  48  for  Pr  =  6.0  shows  a 
definite  increase  in  the  slope  of  the  correlation  r-6  to  —  y~4.  The  other  quantities 
v \r\ u'.8'.uv.u8  do  not  change  with  Pr.  Only  v8  is  affected.  As  Prandtl  number  is 
increased.  v8  shows  a  gradual  increase  of  slope,  indicating  a  smooth  transition  from 
the  third  power  law  to  the  fourth  power  law. 

The  following  tables  show  values  of  Prj  at  the  wall  for  the  present  compu¬ 
tation  and  three  other  models  which  are  qualitatively  similar.  The  listed  wall  value 
of  Prj  for  the  present  computations  is  the  value  at  y~  —  0.4.  This  is  the  closest 
point  which  is  not  grid  dependent,  as  previously  noted. 
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Table  5.3 

Values  of  PrT  at  the  wall 


Pr 

Present 

Computation 

Pr 

Geshev 

0.72 

0.83 

1.00 

1.00 

1.50 

1.10 

2.00 

1.50 

3.00 

1.76 

5.00 

2.30 

6.00 

4.62 

10.0 

3.10 

Pr 

Maksin 

Pr 

Thomas 

0.70 

0.60 

0.72 

1.00 

2.00 

1.30 

1.00 

1.00 

3.00 

1 .50 

2.00 

1.40 

8.00 

1.80 

5.00 

2.40 

The  wall  values  of  PrT  are  tabulated  for  values  of  Pr  somewhat  different  from  the 
values  used  in  the  present  computations  due  to  the  fact  that  the  results  for  some 
models  were  reported  only  for  certain  Pr  and  analytic  equations  were  not  available. 

For  the  present  model. there  is  a  common  intersection  point  near  y~  =  3 
where  PrT  is  nearly  independent  of  the  Pr.  Only  one  other  theory  exhibits  a 
crossover  point  like  this,  and  that  is  the  model  of  Thomas.  However  in  this  case 
the  crossover  is  not  at  a  common  point  and  is  near  y~  =r  20.  In  the  outer  region 
of  the  viscous  sublayer  (y~  >  25).  the  effect  of  Pr  on  PrT  is  relatively  small.  This 
is  anticipated  since  turbulent  transport  in  this  region  is  greater  than  molecular 
transport.  The  trend  of  the  behavior  of  Prj  with  Pr  in  the  region  away  from 
the  wall  is  consistent  with  the  general  rules  of  Reynolds  and  Jischa  (mentioned  in 
chapter  1)  which  are  Prj  >  1  for  Pr  <  I  and.  Prj-  <  1  for  Pr  >  1.  The  theories  of 
Geshev  and  Maksin  do  not  show  this  behavior. 

The  computed  effect  of  a  mean  streamwise  pressure  gradient  on  PrT  for 
air  (Pr  .  .72)  with  an  isothermal  wall  boundary  condition  is  relatively  small  for 


both  adverse  and  favorable  pressure  gradients,  as  as  illustrated  in  figure  42.  These 

can  be  compared  with  the  computational  and  experimental  data  shown  previouslv 

dp 

in  figure  2.  The  effect  of  —  on  PrT  near  the  wall  is  the  same  as  near  the  outer 

di 

edge.  Near  the  outer  edge,  the  trend  of  increasing  Pr p  with  increasing  favorable 
pressure  gradient  deduced  form  experimental  data  (Launder  1978.  Blackwell  1972). 
and  stated  as  a  general  rule  by  Reynolds  and  Jischa,  is  not  seen.  The  opposite 
trend  of  decreasing  Prj  with  increasing  favorable  pressure  gradient  is  observed  and 
this  agrees  with  the  behavior  indicated  by  the  recent  theory  of  Thomas  (1982). 
The  following  two  tables  illustrate  the  trend  of  decreasing  Pr?  with  with  increasing 
favorable  pressure  gradient  for  the  present  model  and  for  Thomas'  model.  The  first 
table  is  for  the  near  wall  region  (y~  =  5).  and  the  second  for  the  outer  region  of  the 
viscous  sublayer  (y*  =  20).  Both  tables  are  for  Pr  =  0.72.  The  experimental  data 
of  Blackwell,  also  included  in  the  tables,  are  seen  to  show  the  opposite  behavior. 


Table  5.4 


Values  of  Pry 

at  y~  -  5 

dp' 

Present 

Thomas 

Blackwell 

di  * 

Computation 

Computation 

Experiment 

0.01 

0.93 

— 

— 

-0.00 

0.9* 

0.85 

1.80 

-0.01 

1.04 

1.00 

1.60 

Table  5.5 


Values  of  Pr  j 

■at  y~  =20 

dp~ 

Present 

Thomas 

Blackwell 

dx~ 

Computation 

Computation 

Experiment 

-0.01 

0.87 

— 

— 

-0.00 

0.99 

0.90 

0.85 

-0.01 

1.04 

1.10 

0.75 
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Some  of  the  boundary  condition  parameters  were  found  to  effect  the  tem¬ 
perature  result?.  For  example,  the  maximum  temperature  intensity  seems  most 
strongly  influenced  by  the  applied  ut  and  vF  intensity  levels.  The  higher  those  in¬ 
tensity  levels,  the  higher  6'^  ...  The  applied  6(  could  be  set  to  zero,  however,  and 
the  0'rnai  is  still  about  the  same.  Similarly,  the  Rvp  and  Ruti  boundary  conditions' 
main  effect  is  on  the  0(y)  profile,  however,  the  variations  due  to  these  boundary 
conditions  are  small. 

Figures  23  for  Pr  -  0.72  and  34(a)  for  Pr  -  6.0  show  the  mean  tem¬ 
perature  profiles  for  the  different  Prandtl  numbers.  As  the  molecular  Pr  number 
is  increased,  the  temperature  profiles  become  imbedded  closer  to  the  wall  as  is 
expected.  Numerical  problems  occur,  however,  at  Pr  >  6.  Gradients  of  mean 
temperature  start  to  oscillate,  causing  positive  gradients  with  a  heated  wall  when 
the  mean  gradients  should  be  negative.  Adjustment  of  the  boundary  condition  pa¬ 
rameters.  cL-.  ou  o.  helps  a  little  but  at  Pr  >  10.  the  mean  temperature  starts 
developing  major  inflectori  points  in  its  profile.  The  inherent  steep  gradient  of  6 
near  the  outer  boundary  becomes  an  increasingly  major  numerical  problem  as  the 
Pr  increases. 

There  is  considerable  uncertainty  about  the  dependence  of  the  outer  edge 
temperature  boundary  conditions  on  Pr.  In  this  study,  no  dependence  of  the  quan¬ 
tities  6[.  (Ruj)t.  and  (/?v.«)f  on  Pr  is  included.  For  the  mean  temperature  6e  in  the 
temperature  boundary  condition.  Kader’s  empirical  equation  (1981)  was  used  to 
formulate  a  Pr  dependence.  Its  applicability  is  from  .001  <  Pr  <  1000  well  within 
the  Pr  range  of  this  study.  Kader's  equation  is 


2-5(2-^) 

O'  Pry~  exp(  -  I)  -  2. 12  ln((  1  -  y  ) - — — jj~)  - \(Pr)  exp(-  1  T)  (5.2a) 
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\(Pr)  =  (3.85Pr3  -  1.3)2  -  2.12  In  Pr 
r  =  Oi  (Pry-)4 

1  *5  Pr^y- 


(5.2c) 


As  regards  the  /?[W  correlation  at  the  outer  edge.  Lawn  (Launder  1978)  has 
suggested  that  this  correlation  is  independent  of  Pr.  and  that  Pr  effects  are  only  to 
be  seen  in  the  O'  values.  The  present  computed  results  do  not  conform  with  this. 
At  the  outer  boundary  (y~  -  40).  Rvfi  decreases  with  increasing  Pr  as  does  Ru(.. 
Temperature  intensity  6'  increases  with  increasing  Pr.  When  a  Prandtl  number 
dependence  in  the  boundary  conditions  was  attempted,  the  quantities  involving 
temperature  would  vary  slightly.  Because  no  experimental  results  were  available  for 
and  R,  h  at  higher  Pr,  it  was  decided  not  to  try  to  construct  such  a  dependence. 

Fulachier  (1984)  states  that  the  mean  temperature  intensity  is  determined 
by  the  mean  turbulent  kinetic  energy.  This  was  shown  to  be  plausible  by  the 
similarity  of  the  spectral  distributions  of  temperature  and  kinetic  energy.  Some 
computational  runs  were  made  with  a  temperature  boundary  condition  proportional 
to  kinetic  energy. 


6 »  =  cq  (5.3) 

It  was  found  that  c  -  1 .0  gave  good  results  for  mean  temperature  and  tempera¬ 
ture  intensity  very  similar  to  the  two  component  model.  With  this  temperature 
boundary  condition,  there  is  no  control  over  (/?,.* i)f  and  {Rue)e.  With  this  model. 
(/?lf.),  =  0.45  and  (Rufi)f  =  0.98.  Results  for  Prandtl  numbers  up  to  6  gave  similar 
results  to  the  two  component  model  for  temperature.  Fulachier  did  not  include 
any  Prandtl  number  dependence  which  should  also  be  included  in  the  expression 
5.3.  This  shows  that  this  may  be  a  boundary  condition  for  temperature  which  may 
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be  physically  correct.  This  boundary  condition  yielded  results  very  close  t’  e  two 
component  model. 

Behavior  of  the  computed  mean  streamwise  velocity  at  the  wall  is  consistent 

with  experiment  when  there  is  an  imposed  mean  streamwise  pressure  gradient  (See 

figures  19. 37(a). 39(a)).  The  mean  temperature  gradient  is  essentially  unchanged. 

(See  figures  23. 37(e). 39(e)).  These  pressure  gradient  runs  are  for  Pr  =  0.72.  With 

a  favorable  pressure  gradient,  the  mean  velocity  gradient  increased  slightly.  With 

— -  -  -  .01.  the  mean  velocitv  gradient  increased  by  4%  and  the  mean  temperature 
dx 

gradient  increased  by  1.0%.  The  adverse  pressure  gradient  had  the  opposite  effect 

with  the  mean  velocitv.  With  ~  =  .01.  the  mean  velocitv  gradient  decreased  bv 

dx 

11%  and  the  mean  temperature  gradient  increased  by  1.0%.  Experiment  shows  that 
u~  and  0~  are  relatively  insensitive  to  pressure  gradient  for  the  region  y~  <  30.  The 
changes  of  u  ~  and  9~  for  the  present  computation  are  small  and  within  experimental 
data  scatter  (Thomas  and  Benton  1982). 

The  values  of  u'  and  O'  for  Pr  =  0.72  are  considerably  higher  than  experi¬ 
mental  measurements  (See  figures  20  and  24).  These  values  are  a  strong  function  of 
the  applied  intensity  of  rf.  The  high  u'  does  not  seem  to  effect  uv  much.  However. 

since  uv  and  -  both  scale  with  u'.  and  their  ratio  appears  in  the  defining  equa- 
dy 

tion  for  Ptj.  some  of  the  problem  with  u'  will  tend  to  be  cancelled  out.  A  similar 

—  dO 

argument  can  be  made  for  the  high  O'  values  since  v 0  and  —  also  occurs  as  a  ratio 

9y 

in  the  equation  for  Fr j-.(Tennekes  and  Lumley  1972) 

Some  qualitative  statements  can  be  made  about  u' .6',  and  p' .  Using  tw'o- 
dimensional  boundary  layer  analysis  els  outlined  by  Cebeci  and  Bradshaw  (1984). 
two  relations  ran  be  derived  with  the  assumptions 

p'  p'  . 

r  <<  u.u  -  0.  —  <<  —  (o.4a) 

P  P 

which  are 


Ou  j it  ~  0,  uf 


(5.4  b) 
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These  two  relations  show  that  both  O'  and  p'  behave  similarly  to  u' .  Figures 
20.24.30(d)  show  that  this  is  qualitatively  true  for  Pr  =  .72. 

Remarks  on  some  of  the  other  turbulent  quantities  ran  be  made.  The  flow 
being  computed  has  a  two-dimensional  main  flow.  That  means  u  -  f{y).  v  =  0.  and 
u'  -  0.  Figures  27(a)  and  33(b)  show  that  f  and  u ■  are  small  but  not  zero.  Using 
the  homogeneity  and  the  symmetry  condition  of  the  flow  in  planes  perpendicular 
to  the  spanwise  z  direction,  the  correlations  involving  the  spanwise  velocity  tr  have 
to  be  zero.  Figures  28(a). 28(b). 28(e)  show  that  the  correlations  utT.mr.  and  wO  for 
Pr  -■  0.72  are  small  but  again  not  exactly  zero.  It  can  be  seen  that  the  higher  order 
correlations  do  not  match  experimental  values  all  that  well.  This  can  be  seen  by 
looking  at  the  plots  for  the  skewnesses  of  u,0.  and  uT  and  the  flatnesses  of  u.O.  and 
ST.  These  are  presented  in  figures  29(a)  to  29(f)  for  Pr  =  0.72. 

§5.4  Effect  on  Heat  Transfer 

The  computed  Pr?  was  inserted  into  the  one  dimensional  mean  temperature 
equation  (equation  1.5).  This  was  done  to  see  what  effects  the  numerically  computed 
variations  with  molecular  Prandtl  number  would  have  on  the  mean  temperature 
profile  and  the  convective  heat  transfer  coefficient  relative  to  the  corresponding 
effects  of  some  reference  model.  For  this  comparison  the  base  model  used  for  PrT 
is  that  of  Kays  and  Crawford  (1980).  Qualitatively  this  model  is  different  from 
the  present  model  in  that  PrT  for  the  Kavs-CrawTord  model  decreases  with  Pr. 
whereas  it  increases  in  the  present  model.  Also,  this  is  an  established  model  for 
Ptj  known  to  give  good  results  for  temperature  profiles  in  fluids  of  different  Prandtl 
numbers.  The  numerically  computed  PrT  compared  to  the  model  of  Kays-Crawford 
gives  somewhat  different  temperature  profiles  as  shown  in  figures  49-52  and  hence 
would  yield  a  corresponding  difference  in  the  convective  heat  transfer  coefficient 
Again.  only  an  approximate  heat  transfer  coefficient  can  be  obtained  using  the 
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values  of  temperature  at  the  outer  boundary  instead  of  the  freestream  temperature. 
The  differences  A/t  in  heat  transfer  are  5.4%  for  Pr  =  .72.8.2%  for  Pr  =  1.5. 
9.4%  for  Pr  -  3.0.  and  10.7%  for  Pr  —  6.0.  These  differences  are  not  very  large 
and  increase  slightly  with  Pr  for  the  given  Pr  range.  It  must  be  noted  that  the 
temperature  differences  even  for  Pr  —  6.0  are  still  within  the  experimental  error 
range  of  the  temperature  measurements  for  Pr  =  .72  of  about  15%  (see  figure  23). 
All  the  temperature  profiles  using  the  computed  Prr  fall  below  the  profiles  of  the 
Kays-Crawford  model. 

To  see  the  effect  on  heat  transfer  of  the  computed  variations  of  PrT  with 
pressure  gradient,  some  similar  computations  were  made  using  equation  1.5.  Refer- 
ing  to  figure  53.  it  can  be  seen  that  the  effect  on  temperature  of  the  computed 
variations  of  Prj  with  streamwise  pressure  gradient  is  negligible.  The  change  in 
the  convective  heat  transfer  coefficient  is  less  than  1.5%. 

Figures  41  and  42  show  that  for  the  cases  of  varying  molecular  Prandtl 
number  and  pressure  gradient  the  turbulent  Prandtl  number  is  essentially  a  constant 
(0.9)  throughout  most  of  the  viscous  sublayer.  The  only  region  where  this  is  not 
true  is  y'  <  3  and  only  at  Pr  >  1.0.  Once  again,  the  one-dimensional  temperature 
analysis  is  employed  to  compare  the  temperature  profiles  using  the  computed  Prr 
with  those  using  a  constant  Prj  of  0.9.  Figures  54  and  55  show  that,  for  Pr  =  0.72 
and  Pr  -  6.0.  the  temperature  profiles  for  variable  and  constant  Prj  are  very 
similar,  with  differences  in  the  approximate  heat  transfer  coefficient  being  less  than 
5%.  This  should  not  be  surprising  since  the  turbulent  Prandtl  number  variations 
occur  at  a  y~  well  within  the  y~  shown  earlier  to  give  significant  variations  in  the 
temperature  profile  (refer  to  chapter  2).  For  Pr  —  6.0,  this  interface  value  was 
found  to  be  at  y~  7.5.  For  higher  Pr.  the  interface  value  is  even  farther  out. 

It  appears  that  the  assumption  of  constant  PrT  =  0.9  across  the  viscous 
sublayer  works  very  well  for  the  Prandtl  number  range  (0.72-6.0).  Launder  (1978) 
states  that  the  assumption  of  constant  turbulent  Prandtl  number  with  Prr  0.9 
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yields  good  heal  transfer  results  for  a  wide  range  of  Pr  from  0.6  to  2000.  Since 
the  model  of  constant  turbulent  Prandtl  number  gives  temperature  profiles  which 
agree  with  the  temperature  profiles  of  the  present  model  well  within  experimental 
uncertainty  (=s  5^7).  it  appears  adequate  in  practical  heat  transfer  computations  to 
neglect  the  variations  in  Ptj  across  the  viscous  sublayer. 

§5.5  Effect  of  a  Second  Component  in  the  ve  Boundary  Condition 

For  one  set  of  computations,  a  second  component  was  added  to  the  te 
boundary  conditions  to  represent  medium  scale  eddies  (MSE).  The  ve  boundary 
condition  was 


r,  -  sinf.V]  -  ovi)Tsinc  -  \  2(J*  -  J?)  sin(3.\'u2T  -  ovj)sin(M  (5.5) 

This  provides  a  more  realistic  value  for  ~  which  experimentally  is  measured  to  be 
about  0.6  (Ota  and  Chapman  1983).  With  only  a  SSE  component  it  is  forced  to  be 
1.0  Additionally  this  second  component  allows  for  more  free  parameters  for  adjust¬ 
ing  the  computational  values.  When  this  model  was  run.  however,  problems  with 


generated  pre--u 


r<  wave-  were  increased  to  the  point  where  the  temperature  profile 


wa-  m  vcn'h  affected  in  an  unrealistic  manner.  The  additional  component  aggra- 
vated  the  boundar>  c ondit ion  problem  at  the  outer  edge.  This  seems  to  confirm  the 
fact  that  the  high  pressure  intensity  is  related  lo  the  oscillating  vf  boundary  con¬ 
dition  Ndditional  MSF  components  were  added  to  the  other  boundary  conditions. 
u..u..  and  b.  to  t r \  to  alleviate  the  severe  pressure  build-up.  These  additional 
component-  did  not  improve  the  temperature  profile,  and  aggravated  the  pressure 
build-up  The  energ\  boundary  condition  with  an  imposed  oscillating  v(  at  the 
boundarv  would  have  to  fie  improved  before  a  more  complex  (realistic)  r£  bound- 
arc  condition  can  be  applied. 

y5.6  Solid-Fluid  Interface  Boundary  Condition 

Thus  far  an  isothermal  wall  boundary  condition  has  been  employed,  in  part 
because  of  its  simplicitv  in  implement  at  ion  and  also  because  most  theories  and 
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experiments  for  turbulent  heat  transfer  assume  this  boundary  condition.  Such  a 
boundary  condition  is  not  exact  in  physical  experiments,  since  wall  temperature 
fluctuations  can  be  measured  (Meek  and  Baer  1973).  Meek  and  Baer  have  con¬ 
cluded  that  that  the  wall  temperature  boundary  condition  is  very  important  in 
determining  turbulent  heat  transfer  physics.  For  this  reason  a  more  complete  spec¬ 
ification  of  the  solid-fluid  temperature  boundary  condition  was  investigated,  and 
used  as  implicit  boundary  conditions  in  the  YSL  computation.  The  conventional 
solid-fluid  interface  boundary  conditions  were  used  to  see  if  the  abnormally  high 
temperature  intensity  peak  wa*-  due  to  the  imposition  of  the  isothermal  wall  tem¬ 
perature  boundary  condition. 

The  solid-fluid  boundary  conditions  for  continuity  of  temperature  and  heat 

flux  are 

0[y~  -  o =  0{y~  =  0 j  (3.6a) 

,  do  do  _ 

(y  0)  /;L.j  -  --  -  (y  =  0)*cltd  (a. 66) 

dy  dy 

With  these  boundary  conditions,  a  stable  algorithm  similar  to  the  present  one  was 
developed  by  Degani  and  Brosh  11983).  This  solid-fluid  interface  boundary  con¬ 
dition  was  run  with  three  different  material  interfaces:  air-copper. air-glass,  and 
water-copper.  The  computed  values  are  reasonable  compared  to  the  experimentally 
measured  values  of  6',,  of  Meek  and  Baer,  however  the  values  cannot  be  compared 
directly  "ince  their  fluid  was  an  oil  (tetralin)  which  has  a  Pr  number  beyond  the 
stable  range  for  the  present  computational  code.  Meek  and  Baer  had  6'u,  ^  0.2  for 
their  oil  flow.  The  following  table  shows  the  values  of  6'u  for  the  different  interfaces. 

Table  5.6 
Values  of  O' 

_  _ _ _ _ U 

Interface 
air  copper 
air  glass 
water  copper 

4* 


0.027 

0.053 

0.112 


It  can  be  seen  in  figures  56(a)-56(c)  that  use  of  the  solid-fluid  interface  condition 
does  not  have  a  major  effect  on  the  temperature  intensity.  Thus,  wall  tempera¬ 
ture  fluctuations  seem  to  have  only  a  very  localized  effect  when  compared  to  the 
isothermal  results.  The  relative!)  high  peak  computed  for  temperature  intensity, 
therefore,  does  not  seem  to  be  a  consequence  of  the  wall  temperature  boundary 
condition.  Other  turbulence  terms  involving  temperature  quantities  are  similarly 
unaffected  in  the  region  away  from  the  wall  in  the  presence  of  wall  temperature  fluc¬ 
tuations.  The  wall  temperature  fluctuations  do  not  seem  to  have  the  importance 
that  Meek  and  Baer  have  suggested.  It  is  believed  that  the  high  temperature  inten¬ 
sity  is  due  to  some  unrealistic  aspect  of  the  outer  boundary  condition  at  y  =40 
and  not  due  to  the  isothermal  wall  temperature  boundary  condition. 
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6.  Concluding  Remarks 

\  tun  component  model  has  been  developed  for  computing  the  distribution 

of  t  urbulen*  Prandt  I  number  in  the  compressible  viscous  sublayer  for  flows  of  varying 

mean  streamuiM-  pressure  gradient  and  varying  molecular  Prandtl  number.  The 

principal  problem  encountered  involved  finding  a  set  of  boundary  conditions  that 

yields  stable  numerical  solutions  without  numerical  decoupling.  A  comparison  of 

computed  results  with  available  experiments  has  been  made  for  the  four  quantities 
du  dB. 

(uv.va.  -~- .  — )  that  determine  PrT.  The  agreement  is  satisfactory,  although  it  is 
oy  oy 

noted  that  data  presently  available  for  v6  exist  only  near  the  outer  edge  of  the 
viscous  sublayer. 

The  range  of  values  for  the  computed  Ptt  is  about  in  the  middle  of  the  wide 

range  of  predicted  values  of  Pr ?  from  various  experiments  and  existing  theories. 

The  theories  of  Thomas,  Maksin.and  Geshev  are  qualitatively  consistent  with  the 

present  computed  trends  of  the  variations  of  Pr t  with  molecular  Pr.  All  the  rest 

surveyed  (6  theories)  show  differing  trends. 

There  is  a  strong  dependence  of  Prj  on  Pr  near  the  wall  (y+  <  3).  but 

the  dependence  of  Prj  on  —  is  verv  slight  throughout  the  viscous  sublayer.  For 

ox 

y~  >  3.  the  turbulent  Prandtl  number  is  essentially  a  constant,  in  the  range  0. 9-1.0. 
Th-  small  effect  of  mean  pressure  gradient  on  Ptj  is  consistent  w;ith  the  available 
experimental  data.  The  behavior  of  computed  Prj  with  pressure  gradient  agrees 
with  the  trend  predicted  by  Thomas.  The  present  computation  agrees  qualitatively 
with  the  computational  model  of  Thomas  both  for  the  case  of  varying  Prandtl 
number  and  pressure  gradient. 

Molecular  Prandtl  numbers  between  0.7  and  6.0  affect  the  turbulent  Prandtl 
number  only  in  a  region  so  close  to  the  wall  that  the  temperature  profiles  and 
heat  transfer  are  essentially  unaffected.  In  view  of  this,  and  of  the  small  effect 
of  streamwise  pressure  gradient  on  turbulent  Prandtl  number,  the  assumption  of 
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constant  Ptj  =  0.9  across  the  viscous  sublayer  appears  adequate  for  practical 
computations  of  heat  transfer. 

It  has  been  shown  that  the  temperature  profile  and  convective  heat  transfer 
coefficient  computed  using  the  present  Pr 7-  model  for  both  varying  mean  stream- 
wise  pressure  gradient  and  molecular  Prandtl  number  differ  only  slightly  from  the 
temperature  profiles  computed  with  the  Pr?  model  of  Kays  and  Crawford.  The 
small  temperature-profile  differences  between  the  two  models  are  somewhat  larger 
for  the  cases  of  varying  Pr  than  for  varying  pressure  gradient. 

A  solid-fluid  interface  boundary  condition  was  investigated  in  addition  to 
an  isothermal  wall  boundary  condition  to  see  if  the  wall  temperature  boundary 
condition  was  critical  in  determining  the  fluctuating  temperature  intensity.  It  was 
found  the  the  solid-fluid  interface  condition  produced  wall  temperature  fluctuations 
as  expected,  however,  their  effect  was  not  large  and  was  very  localized  near  the 
wall.  Therefore,  it  is  concluded  that  the  wall  temperature  boundary  condition  is 
not  a  critical  factor  in  determining  the  intensity  of  temperature  fluctuations  across 
the  viscous  sublayer. 

The  main  problem  encountered  in  developing  the  computational  model 
stems  from  the  numerical  boundary  condition  for  energy  in  the  presence  of  an 
oscillating  normal  velocity  at  the  outer  edge  boundary.  If  this  problem  can  be  fully 
resolved,  a  more  versatile  model  could  be  developed  using  the  present  compressible- 
flow  algorithm.  The  numerical  problems  encountered  in  this  research  are  believed 
to  be  representative  of  those  to  be  encountered  when  a  compressible  code  is  run  at 
low  Mach  number  and  a  central  differencing  algorithm  is  run  without  dissipation. 
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APPENDIX  I 

Decomposition  of  Jacobian  and  Viscous  Matrices 
The  Jacobian  and  viscous  matrices  (B  -  J~]BXJ)  =  dF  dq  -  dFx  dq  and 

(C  -  J~1CXJ)  =  dG  dq  -  dGv  dq  obtained  in  the  time  linearization  of  F  -  Fx 
and  G  -  Gv  are  each  decomposed  into  two  matrices  .  one  (subcript  0)  that  does  not 
contain  any  cross  derivative  terms  and  another(subcript  x)  that  contains  all  such 


terms. 
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APPENDIX  II 

Temperature  Solution  for 
Oscillating  Shear  Flow 

An  exact  analytical  solution  has  been  derived  for  oscillating  shear  flow 
(OSF)  with  heat  transfer  in  an  incompressible  fluid  with  constant  properties.  Since 
the  momentum  equations  are  independent  of  the  thermal  energy  equation  in  incom¬ 
pressible  flow,  the  velocity  field  is  the  same  as  for  the  case  without  heat  transfer 
derived  by  Chapman  and  Kuhn  (1981).  With  y  as  the  coordinate  normal  to  the 
infinite  flat  surface,  and  u  and  r  as  the  velocity  components  in  the  streamwise  and 
spanwise  directions,  respectively,  we  have 

u  ~  .4  sinnf  -  e~k"ys\n {nt  -  kny)\  -  cy 
iv  ~  B  sin(mt  -  o)  -  e~krny s\n(mt  -  kmy) 


where 


,  n  .  _  ;  ™ 

"  "  \  2i  m~  V  2 * 


The  oscillating  frequencies  are  n  in  the  streamwise  direction  and  m  in  the  spanwise 
direction,  and  the  mean  velocity  profile  is  u  -  cy.  For  this  velocity  field,  the  solution 
to  the  thermal  energy  equation 
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Molecular  Pr  -  .72  (Air) 


perimental  represent  all 
ith  an  adverse  pressure 


Figure  3  Theoretical  uncertainty  envelope  for  the  turbulent  Prandtl  number  for  Pr 
(Air)  and  zero  pressure  gradient.  Eight  models  for  Pr-p  shown. 
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Theoretical  uncertainty  envelope  for  the  turbulent  Prandtl  numb* 
(Water)  and  zero  pressure  gradient.  Seven  models  for  PrT  shown 
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Figure  6  Theoretical  uncertainty  envelope  for  the  turbulent  Prandtl  number  for  various 
molecular  Prandtl  number  and  zero  pressure  gradient.  Two  models  for  PrT 
shown. 
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Figure  7  Temperature  profiles  for  a  two-layer  step  model  of  the  turbulent  Prandtl  number 
and  zero  pressure  gradient.  This  shows  the  effect  of  the  value  of  Prj  near  the 
wall.  For  Pr  —  0.72,  value  of  Prj  is  switched  at  y+  =  13.7.  Ah  is  an  approximate 
change  in  the  convective  heat  transfer  coefficient  in  per  cent. 
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Figure  8  Temperature  profiles  for  a  two-layer  step  model  of  the  turbulent  Prandtl  number 
and  zero  pressure  gradient.  This  shows  the  effect  of  the  value  of  Prr  near  the 
wall.  For  Pr  —  6.0,  value  of  Prj  is  switched  at  y+  =  7.5.  Ah  is  an  approximate 
change  in  the  convective  heat  transfer  coefficient  in  per  cent. 


73 


LSE 


Figure  12  Spanwise  view  of  computational  domain  for  VSL  model  with  contra-rotating  vor 
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'ode  verification  test  with  oscillating  shear  flow  and  an  isothermal  wall  tern- 


Figure  16(c)  Contour  plot  of  the  vertical  velocity  component  v  before  any  averaging  is  applied 
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Figure  18  Grid  dependency  illustrated  near  wall  with  Prj  as  grid  is  refined. 
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Figure  J9  Computed  mean  streamwise  velocity  d*  for  Pr  =  0.72, isothermal  wall  tempera 
ture  boundary  condition, zero  mean  streamwise  pressure  gradient.  Symbols  rep 
resent  experimental  data. 
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Figure  21  Computed  Reynolds  stress  uT~  for  Pr  =  0.72, isothermal  wall  temperature  bound 
ary  condition. zero  mean  streamwise  pressure  gradient.  Symbols  represent  exper 
imental  data. 
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Figure  22  Computed  correlation  coefficient  Ruv  for  Pr  =  0.72, isothermal  wall  temperature 
boundary  condition, zero  mean  streamwise  pressure  gradient.  Symbols  represent 
experimental  data. 
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Figure  23  Computed  mean  temperature  for  Pr  -  0.72, isothermal  wall  temperature 
boundary  condition. zero  mean  streamwise  pressure  gradient.  Symbols  represent 
experimental  data. 
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Figure  25  Computed  correlation  coefficient  Rue  for  Pr  =  0.72, isothermal  wall  temperature 
boundary  condition, zero  mean  streamwise  pressure  gradient.  Symbols  represent 
experimental  data. 
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Figure  26  Computed  correlation  coefficient  Rve  for  Pr  —  0.72, isothermal  wall  temperature 
boundary  condition, zero  mean  streamwise  pressure  gradient.  Symbols  represent 
experimental  data. 


Figure  27  (a)  Mean  normal  velocity  v,(b)  Mean  spanwise  velocity  tt\(c)  Normal  velocin 
intensity  r'.o  =  Schildnecht  (1979),?  =  Laufer  (1953),  (d)  Spanwise  velocity  in¬ 
tensity  w'.  Symbols  same  as  (c),  for  Pr  =  0.72  and  zero  pressure  gradient. 


igure  28  (a)  Correlation  uu\(b)  Correlation  ru\(c)  Correlation  u 6.  (d)  Correlation  v6 ,  (e) 
Correlation  tt'0,(f)  Correlation  p$,  for  Pr  =  0,72  and  zero  pressure  gradient. 
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Figure  29  (a)  Skewness  of  u,  o  =  Ueda  (1975),  o  =  Elena  (1979),  (b)  Flatness  of  u,  Symbols 
same  as  (a),  (cj  Skewness  of  6,  Symbols  same  as  (a),  (d)  Flatness  of  6,  Symbols 
same  as  (a),(e)  Skewness  of  uu,  x  =  Gupta  (1972),  (f)  Flatness  of  uv,  Symbol 
same  as  (e),  for  Pr  =  0.72  and  zero  pressure  gradient. 
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Figure  32  (a)  Correlation  il-6,  (b)  Correlation  pd.  (c)  Mean  pressure  p,  (d)  Pressure  intensity 
p' .  (e)  Mean  density  p,  (f)  Density  intensity  p\  (g)  Skewness  of  temperature  5. 
(h)  Flatness  of  temperature  FV,  for  Pr  =  1.5  and  zero  pressure  gradient. 
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Figure  33  (a)  Mean  temperature  0.(b)  Temperature  intensity  ff\(c)  Correlation  coefficient 
Ru*.  (d)  Correlation  coefficient  R (e)  Correlation  uO,  (f)  Correlation  16.  for 
Pr  —  3.0  and  zero  pressure  gradient. 
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Figure  34  (a)  Correlation  uff,  (b)  Correlation  p9.  (c)  Mean  pressure  p.  (d)  Pressure  intensity 
p',  (e)  Mean  density  p.  (f)  Density  intensity  p',  (g)  Skewness  of  temperature  S. 
(h)  Flatness  of  temperature  F$ ,  for  Pr  ~  3.0  and  zero  pressure  gradient. 
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Figure  37  (a)  Mean  streamwise  velocity  u,  (b)  Streamwise  velocity  intensity  u  ,  (c)  Reynolds 
stress  -uT,  (d)  Reynolds  stress  correlation  coefficient  -Ru i.  (e)  Mean  temper¬ 
ature  0.(f)  Temperature  intensity  0',(g)  Correlation  coefficient  -Rue-  (h)  Corre- 

dp 

lation  coefficient  Rve ,  for  Pr  =  0.72  and  —  —  .01,  o  =  zero  pressure  gradient 

ox 

case. 
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(a)  Mean  streamwise  velocity  u.  (b)  Streamwise  velocity  intensity  u  ,  (c)  Reynolds 
stress  -  tub  (d)  Reynolds  stress  correlation  coefficient  ~/?uv,  (e)  Mean  tempera¬ 
ture  0.(f)  Temperature  intensity  0',(g)  Correlation  coefficient  - Ru# ,  (h)  Correla¬ 
te 

tion  coefficient  Rtt>.  for  Pr  =  0.72  and  —  =  -.01,  c  =  zero  pressure  gradient 
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Figure  49  Temperature  profile  comparison  between  the  Kays-Crawford  model  of  Prj  and 
the  computed  profile  of  Prj-  for  Pr  =  0.72.  A/i  represents  approximate  change 
in  the  convective  heat  transfer  coefficient. 
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Figure  51  Temperature  profile  comparison  between  the  .Kays-Crawford  model  of  PrT  and 
the  computed  profile  of  Prj  for  Pr  -  3.0.  Ah  represents  approximate  change  in 
the  convective  heat  transfer  coefficient. 
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Figure  52  Temperature  profile  comparison  between  the  Kays-Crawford  model  of  PrT  and 
the  computed  profile  of  Prj-  for  Pr  =  6.0.  A/i  represents  approximate  change  in 
the  convective  heat  transfer  coefficient. 
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Figure  54  Temperature  profile  comparison  between  a  constant  Prj  =  0.9  and  the  computed 
profile  of  Prj  for  Pr  =  0.72 represents  approximate  change  in  the  convective 
heat  transfer  coefficient. 
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Figure  55  Temperature  profile  comparison  between  a  constant  Pr t  =  0.9  and  the  computed 
profile  of  PrT  for  Pr  =  6.0.  .  A/i  represents  approximate  change  in  the  convective 
heat  transfer  coefficient. 


